Parallel Variable Distribution Algorithm for Constrained Optimization with Nonmonotone Technique

Abstract

A modified parallel variable distribution (PVD) algorithm for solving large-scale constrained optimization problems is developed, which modifies quadratic subproblem 𝑄𝑃_𝑙 at each iteration instead of the 𝑄𝑃⁰_𝑙 of the SQP-type PVD algorithm proposed by C. A. Sagastizábal and M. V. Solodov in 2002. The algorithm can circumvent the difficulties associated with the possible inconsistency of 𝑄𝑃⁰_𝑙 subproblem of the original SQP method. Moreover, we introduce a nonmonotone technique instead of the penalty function to carry out the line search procedure with more flexibly. Under appropriate conditions, the global convergence of the method is established. In the final part, parallel numerical experiments are implemented on CUDA based on GPU (Graphics Processing unit).

1. Introduction

In this paper, we consider the following nonlinear programming problem: (𝑥),(𝑥)min𝑓s.t.𝑐≤0, (1) where 𝑓∶𝑅^𝑛→𝑅 and 𝑐(𝑥)∶𝑅^𝑛→𝑅^𝑚 are all continuously differentiable. Suppose the feasible set 𝑋 in (1) is described by a system of inequality constraints: 𝑍=𝑥∈𝑅^𝑛∣𝑐_𝑖(𝑥)≤0,𝐼={𝑖∣𝑖=1,2,…,𝑚}. (2)

In this paper, we give a new algorithm based on the method in [1], which partitions the problem variables 𝑥∈𝑅^𝑛 into 𝑝 blocks 𝑥₁,…,𝑥_𝑝, such that 𝑥=𝑥₁,…,𝑥_𝑝,𝑥_𝑙∈𝑅^𝑛𝑙,𝑙=1,…,𝑝,^𝑝_𝑙=1𝑛_𝑙=𝑛. (3)

And assume that 𝑋 has block-separable structure. Specifically, 𝑐(𝑥)=𝑐₁𝑥₁,𝑐₂𝑥₂,…,𝑐_𝑝𝑥_𝑝,𝑐_𝑙𝑥_𝑙∶𝑅^𝑛𝑙⟶𝑅^𝑚𝑙,𝑙∈1,2,…,𝑝,^𝑝_𝑙=1𝑚_𝑙=𝑚. (4)

Parallel variable distribution (PVD) algorithm for solving optimization problems was first proposed by Ferris and Mangasarian [2], in 1994, in which the variables are distributed among 𝑝 processors. Each processor has the primary responsibility for updating its block of variables while allowing the remaining secondary variables to change in a restricted fashion along some easily computable directions. The distinctive novel feature of this algorithm is the presence of the “forget-me-not” term which allows for a change in “secondary” variables. This makes PVD fundamentally different from the block Jacobi [3] and coordinate descent [4] methods. The forget-me-not approach improves robustness and accelerates convergence of the algorithm and is the key to its success. In 1997, Solodov [5] proposed useful generalizations that consist, for the general unconstrained case, of replacing exact global solution of the subproblems by a certain natural sufficient descent condition, and, for the convex case, of inexact subproblem solution in the PVD algorithm. Solodov [6] proposed an algorithm applied the PVD approach to problems with general convex constraints directly use projected gradient residual function as the direction and show that the algorithm converges, provided that certain conditions are imposed on the change of secondary variables. In 2002, Sagastizábal and Solodov [1] proposed two new variants of PVD for the constrained case. Without assuming convexity of constraints, but assuming block-separable structure, they showed that PVD subproblems can be solved inexactly by solving their quadratic programming approximations. This extends PVD to nonconvex (separable) feasible sets and provided a constructive practical way of solving the parallel subproblems. For inseparable constraints, but assuming convexity, they developed a PVD method based on suitable approximate projected gradient directions. The approximation criterion was based on a certain error bound result, and it was readily implementable. In 2011, Zheng et al. [7] gave a parallel SSLE algorithm, in which the PVD subproblems are solved inexactly by serial sequential linear equations, for solving large-scale constrained optimization with block-separable structure. Without assuming the convexity of constraints, the algorithm is proved to be globally convergent to a KKT point. Han et al. [8] proposed an asynchronous PVT algorithm for solving large-scale linearly constrained convex minimization problems with the idea of [9] in 2009, which based on the idea that a constrained optimization problem is equivalent to a differentiable unconstrained optimization problem by introducing the Fischer Function. And in particular, different from [9] the linear rate of convergence does not depend on the number of processors.

In this paper, we use [1] as our main reference on SQP-type PVD method for problems (1) and (4). Firstly, we introduce the algorithm in [1] simply. The original problem is distributed into 𝑝 parallel subproblems which treatment among 𝑝 parallel processors. The algorithm of [1] may result in that the linear constraints in quadratic programming subproblems are inconsistent or the solution of quadratic programming subproblems is unbounded, so that the algorithm may fail. This drawback has been overcome by many researchers, such as [10]. In [1], exact penalty function is used as merit function to carry out the line-search. For the penalty function, as pointed out by Fletcher and Leyffer [11], the biggest drawback is that the penalty parameter estimates could be problematic to obtain. To overcome this drawback, we use a nonmonotone technique to carry out the line search instead of penalty function.

Recent research [12–14] indicates that the monotone line search technique may have some drawbacks. In particular, enforcing monotonicity may considerably reduce the rate of convergence when the iteration is trapped near a narrow curved valley, which can result in very short steps or zigzagging. Therefore, it might be advantageous to allow the iterative sequence to occasionally generate points with nonmonotone objective values. Grippo et al. [12] generalized the Armijo rule and proposed a nonmonotone line search technique for Newton’s method which permits increase in function value, while retaining global convergence of the minimization algorithm. In [15], Sun et al. give several nonmonotone line search techniques, such as nonmonotone Armijo rule, nonmonotone Wolfe rule, nonmonotone F-rule, and so on. Several numerical tests show that the nonmonotone line search technique for unconstrained optimization and constrained optimization is efficient and competitive [12–14, 16]. Recently, [17] gives a method to overcome the drawback of zigzagging with nonmonotone line search technique to determine the step length, which makes the algorithm more flexible.

In this paper, we combine [1] with the ideas of [10, 17] and propose an infeasible SQP-type Parallel Variable Distribution algorithm for constrained optimization with nonmonotone technique.

The paper is organized as follows. The algorithm is presented in Section 2. In Section 3, under mild assumptions, some global convergence results are proved. The numerical results are shown in Section 4. And conclusions are given in the last section.

2. A Modified SQP-Type PVD Method with Nonmonotone Technique

To describe our algorithm, we first give the modified quadratic subproblems of SQP method. Given 𝑥^𝑘∈𝑅^𝑛, we distribute it into 𝑝 blocks, as the iterate 𝑥^𝑘_𝑙∈𝑅^𝑛𝑙,𝑙=1…,𝑝,𝐻_𝑙∈𝑅^{𝑛𝑙×𝑛𝑙},𝑙=1,2,…,𝑝 denote an approximate Hessian.

We use 𝑧_𝑙∈𝑅 to perturb the subproblem of 𝑄𝑃⁰_𝑙 of [1] and give a new subproblem 𝑄𝑃_𝑙 instead of it. Consider 𝑄𝑃_𝑙𝑥^𝑘_𝑙,𝐻^𝑘_𝑙=min_{(𝑑𝑙,𝑧𝑙)∈𝑅𝑛𝑙+1}𝑧_𝑙+12𝑑^𝑇_l𝐻^𝑘_𝑙𝑑_𝑙s.t.𝑔^𝑘_𝑙^𝑇𝑑_𝑙≤𝑧_𝑙,𝑐_𝑖𝑥^𝑘_𝑙+∇𝑐_𝑖𝑥^𝑘_𝑙^𝑇𝑑_𝑙≤𝑧_𝑙,𝑖∈𝐼_𝑙=1,…,𝑚_𝑙. (5)

For convenience, given 𝑥, the 𝑙th block of ∇𝑓(𝑥) will be denoted by 𝑔_𝑙=𝐼_{𝑛𝑙×𝑛𝑙}0_{𝑛𝑙×𝑛𝑙}∇𝑓(𝑥). (6)

In (5), 𝑔^𝑘_𝑙 is 𝑔_𝑙 in step 𝑘. Note that (5) is always feasible for 𝑑_𝑙=0,𝑧_𝑙=max_{𝑖∈𝐼𝑙}{𝑐_𝑖(𝑥^𝑘_𝑙),0}. To test whether constraints are satisfied or not, we denote the violation function ℎ(𝑥) as follows: ℎ(𝑥)=𝑐(𝑥)⁺, (7)where 𝑐_𝑖(𝑥)⁺=max{𝑐_𝑖(𝑥),0}, 𝑖∈𝐼, ‖⋅‖ denotes the Euclidean norm on 𝑅^𝑚. It is easy to see that ℎ(𝑥)=0 if and only if 𝑥 is a feasible point (ℎ(𝑥)>0 if and only if 𝑥 is infeasible).

We describe our modified PVD algorithm as follows.

Algorithm A. Consider the following.

Step 0. Start with any 𝑥₀∈𝑅^𝑛. Choose parameters 𝑢∈(0,1/2), 𝛽∈(1/2,1), 𝛾∈(0,1), positive integer 𝑀>0, and positive definite matrices 𝐻⁰_𝑙∈𝑅^{𝑛𝑙×𝑛𝑙}, 𝑙∈{1,2,…,𝑝}. Set 𝑘=0, 𝑚(𝑘)=0.

Having 𝑥_𝑘, check a stopping criterion. If it is not satisfied, compute 𝑥_𝑘+1 as follows.

Step 1. Parallelization. For each processor 𝑙∈{1,2,…,𝑝}, solve subproblem (5) to get (𝑑^𝑘_𝑙,𝑧^𝑘_𝑙). If 𝑑^𝑘_𝑙=0, 𝑙=1,…,𝑝 stop, otherwise return to Step 2.

Step 2. Synchronization. Set 𝑑_𝑘=𝑑^𝑘₁,…,𝑑^𝑘_𝑝,𝑔^𝑇_𝑘=𝑔^𝑘₁^𝑇,…,𝑔^𝑘_𝑝^𝑇,𝐻_𝑘=𝐻^𝑘₁⋱𝐻^𝑘_𝑝. (8)
If 𝑔^𝑇_𝑘𝑑_𝑘≤−(1/2)𝑑^𝑇_𝑘𝐻_𝑘𝑑_𝑘, return to Step 3; otherwise, set 𝜆_𝑘=1 and return to Step 4.

Step 3. Choose 𝜆_𝑘 which is the largest one in the sequence {1,𝛾,𝛾²,…} satisfying: 𝑓𝑥_𝑘+𝜆𝑑_𝑘≤max_{(𝑘)0≤𝑗≤𝑚}𝑓𝑥_𝑘−𝑗−𝜆𝜇𝑑^𝑇_𝑘𝐻_𝑘𝑑_𝑘. (9)

Step 4. Set 𝑥_𝑘+1=𝑥_𝑘+𝜆_𝑘𝑑_𝑘. Let ℎ=max_{0≤𝑗≤𝑚(𝑘)}[ℎ(𝑥_𝑘−𝑗)]. If ℎ(𝑥_𝑘+1)≤𝛽ℎ, then set 𝑚(𝑘+1)=min{𝑚(𝑘)+1,𝑀}, update 𝐻^𝑘_𝑙 to 𝐻_𝑙^𝑘+1,𝑙=1,…,𝑝, set 𝑘=𝑘+1, and go to Step 1. Otherwise, let 𝑘=𝑘+1, call Restoration Algorithm of [17] to obtain 𝑥^𝑟_𝑘, let 𝑥_𝑘=𝑥^𝑟_𝑘, 𝑚(𝑘)=𝑚(𝑖), and go to Step 1.

Remark 1. In Step 3 of Algorithm A, the nonmonotone parameter 𝑚(𝑘) satisfies 𝑚(0)=0,0≤𝑚(𝑘)≤min{𝑚(𝑘−1)+1,𝑀},𝑘≥1. (10)
For the convenience, we denote 𝑓(𝑥_𝑙(𝑘))=max_{0≤𝑗≤𝑚(𝑘)}[𝑓(𝑥_𝑘−𝑗)], where 𝑘−𝑚(𝑘)≤𝑙(𝑘)≤k.

In a restoration algorithm, we aim to decrease the value of ℎ(𝑥) more precisely, we will use a trust region type method to obtain ℎ(𝑥_𝑘)→0, 𝑘→∞ by the help of the nonmonotone technique. Let 𝑀^𝑖_𝑘(𝑑)=ℎ𝑥^𝑖_𝑘−ℎ𝑥^𝑖_𝑘,Ψ+𝑑_𝑖𝑘(𝑑)=ℎ𝑥^𝑖_𝑘−𝑐𝑥^𝑖_𝑘+∇𝑐𝑥^𝑖_𝑘^𝑇𝑑⁺, (11)

and 𝑟^𝑖_𝑘=𝑀^𝑖_𝑘(𝑑)/Ψ^𝑖_𝑘(𝑑).

Algorithm B is similar to the restoration phase given by Su and Yu [17], we describe the Restoration Algorithm as follows.

Algorithm B. Consider the following.

Step 0. Assume 𝑥⁰_𝑘=𝑥_𝑘, Δ⁰>0, 𝑖=0, 𝜂∈(0,1), 𝑚(𝑖)=0.

Step 1. If ℎ(𝑥^𝑖_𝑘)≤𝛽ℎ, then 𝑥^𝑟_𝑘=𝑥^𝑖_𝑘 stop.

Step 2. Compute maxΨ^𝑖_𝑘(𝑑)s.t.‖𝑑‖≤Δ^𝑖 (12)
to get 𝑑^𝑖_𝑘. Calculate 𝑟^𝑖_𝑘.

Step 3. If 𝑟^𝑖_𝑘≤𝜂, then let 𝑥_𝑘^𝑖+1=𝑥^𝑖_𝑘, Δ^𝑖+1=Δ^𝑖/2, 𝑖=𝑖+1, 𝑚(𝑖)=min{𝑚(𝑖−1)+1,𝑀} and go to Step 2.

Step 4. If 𝑟^𝑖_𝑘>𝜂, then let 𝑥_𝑘^𝑖+1=𝑥^𝑖_𝑘+𝑑^𝑖_𝑘, Δ^𝑖+1=2Δ^𝑖, 𝑖=𝑖+1, 𝑚(𝑖)=min{𝑚(𝑖−1)+1,𝑀} and go to Step 1.

3. The Convergence Properties

To prove the global convergence of Algorithm A, we make the following assumptions.

Assumptions

(A1) The iterate {𝑥_𝑘} remains in compacted subset 𝑆⊂𝑅^𝑛.(A2) The objective function 𝑓 and the constraint functions 𝑐_𝑗(𝑗∈𝐼) are twice continuously differentiable on 𝑅^𝑛.(A3) For all 𝑥_𝑙∈𝑅^𝑛𝑙, the set {∇𝑐_𝑖(𝑥_𝑙)∶𝑖∈𝐼_𝑙(𝑥_𝑙)} is linearly independent, where 𝐼_𝑙(𝑥_𝑙)={𝑖∈𝐼_𝑙∶𝑐_𝑖(𝑥_𝑙)=Φ_𝑙(𝑥_𝑙),Φ_𝑙(𝑥_𝑙)=max_{𝑖∈𝐼𝑙}{𝑐_𝑖(𝑥_𝑙),0}}.(A4) There exist two constants 0<𝑎≤𝑏 such that 𝑎‖𝑑_𝑙‖²≤𝑑^𝑇_𝑙𝐻^𝑘_𝑙𝑑_𝑙≤𝑏‖𝑑_𝑙‖², 𝑙=1,…,𝑝 for all iteration 𝑘 and 𝑑_𝑙∈𝑅^𝑛𝑙.(A5)The solution of problem (12) satisfies Ψ^𝑖_𝑘(𝑑)=ℎ𝑥^𝑖_𝑘−𝑐𝑥^𝑖_𝑘+∇𝑐𝑥^𝑖_𝑘^𝑇𝑑⁺≥𝛽₂Δ^𝑖minℎ𝑥^𝑖_𝑘,Δ^𝑖, (13)where 𝛽₂ is a constant.

Remark 2. Assumptions (A1) and (A2) are the standard assumptions. (A3) is the LICQ constraint qualification. (A4) plays an important role in obtaining the convergence results. (A5) is the sufficient reduction condition which guarantees the global convergence in a trust region method. Under the assumptions, 𝑓 is bounded below and the gradient function 𝑔_𝑙,𝑙=1,…,𝑝 is uniformly continuous in 𝑆_𝑙, where 𝑆=(𝑆₁,…,𝑆_𝑝).

Remark 3. We can use quasi-Newton BFGS methods to update 𝐻^𝑘_𝑙, different from [18], we can use a small modification of 𝑦^𝑘_𝑙 to make 𝐻^𝑘_𝑙 reserve positive definite according to the formula (33) of [17].

Similar to Lemma  1 in [19], we can get the following conclusions.

Lemma 4. Suppose (A1)–(A4) hold, 𝐻^𝑘_𝑙,𝑙=1,…,𝑝 is a symmetric positive definite, then 𝑑^𝑘_𝑙 is well defined and (𝑑^𝑘_𝑙,𝑧^𝑘_𝑙) is the unique KKT point of (5). Furthermore, 𝑑^𝑘_l is bounded over compact subsets of 𝑆_𝑙×𝑃_𝑙, where 𝑃_𝑙 is the set of symmetric positive definite 𝑛_𝑙×𝑛_𝑙 matrices.

Proof. First note that the feasible set for (5) is nonempty, since (𝑑_𝑙,𝑧_𝑙)=(0,max_{𝑖∈𝐼𝑙}{𝑐_𝑖(𝑥_𝑙),0}) is always feasible. It is clear that, if (𝑑^𝑘_𝑙,𝑧^𝑘_𝑙) is a solution to (5) if an only if 𝑑^𝑘_𝑙 solves the unconstrained problem min_{𝑑𝑙∈𝑅𝑛𝑙}12𝑑^𝑇_𝑙𝐻^𝑘_𝑙𝑑_𝑙+max𝑔^𝑘_𝑙^𝑇𝑑_𝑙,max_{𝑖∈𝐼𝑙}𝑐_𝑖𝑥^𝑘_𝑙+∇𝑐_𝑖𝑥^𝑘_𝑙^𝑇𝑑_𝑙,𝑧_𝑙=max𝑔^𝑘_𝑙^𝑇𝑑_𝑙,max_{𝑖∈𝐼𝑙}𝑐_𝑖𝑥^𝑘_𝑙+∇𝑐_𝑖𝑥^𝑘_𝑙^𝑇𝑑_𝑙. (14)(15) Since the function being minimized in (14) is strictly convex and radially unbounded, it follows that 𝑑^𝑘_𝑙 is well defined and unique as a global minimizer for the convex problem (5) and thus unique as a KKT point for that problem. So we have 𝐻^𝑘_𝑙𝑑^𝑘_𝑙+𝑣𝑔^𝑘_𝑙+∇𝑐_𝑙𝑥^𝑘_𝑙^𝑇𝜇^𝑘_𝑙=0,1−𝑣−^𝑚𝑙_𝑖=1𝑢_𝑖=0, (16)(17) due to (17) and 𝑢_𝑖≥0,𝑣≥0, we have 𝑢_𝑖,𝑖∈𝐼_𝑙 is bounded. Since 𝑓,𝑐_𝑖 is twice continuously differentiable and assumption (A1) holds, for all 𝑥_𝑙∈𝑆_𝑙, ‖−(𝑣𝑔_𝑙+∇𝑐_𝑙(𝑥^𝑘_𝑙)^𝑇𝜇^𝑘_𝑙‖ is bounded, set the maximum is 𝑀. For (16), ‖𝐻^𝑘_𝑙𝑑^𝑘_𝑙‖≤𝑀 and combining with Assumption (A4), we have 𝑎‖𝑑^𝑘_𝑙‖²≤𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙≤‖𝑑^𝑘_𝑙‖‖𝐻^𝑘_𝑙𝑑^𝑘_𝑙‖≤𝑀‖𝑑^𝑘_𝑙‖, so that ‖𝑑^𝑘_𝑙‖≤𝑀/𝑎; therefore, 𝑑^𝑘_𝑙 is bounded on 𝑆_𝑙×𝑃_𝑙.

Lemma 5. Suppose that 𝑥^𝑘_𝑙∈𝑅^𝑛𝑙,𝑙=1,…,𝑝,𝐻^𝑘_𝑙 is positive definite matrix, and (𝑑^𝑘_𝑙,𝑧^𝑘_𝑙) is the solution of 𝑄𝑃_𝑙. Then we have the results (B1) and (B2).(B1) The following inequality holds: 𝑧^𝑘_𝑙≤Φ_𝑙𝑥^𝑘_𝑙−12𝑑^𝑘_𝑙^𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙,𝑙=1,…,𝑝, (18) where Φ_𝑙(𝑥_𝑙)=max_{𝑖∈𝐼𝑙}{𝑐_𝑖(𝑥_𝑙),0}.(B2) If 𝑑^𝑘_𝑙=0 then 𝑧^𝑘_𝑙=0,𝑙=1,…,𝑝, and 𝑥_𝑘=(𝑥^𝑘₁,…,𝑥^𝑘_𝑝) is a KKT point of problem (1).

Proof. (B1) Since 𝑑_𝑙=0, 𝑧_𝑙=max_{𝑖∈𝐼𝑙}{𝑐_𝑖(𝑥_𝑙),0} is a feasible point of 𝑄𝑃_l, from the optimality of (𝑑^𝑘_𝑙,𝑧^𝑘_𝑙), we have 𝑧^𝑘_𝑙+12𝑑^𝑘_𝑙^𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙≤𝑧_𝑙+0=max_{𝑖∈𝐼𝑙}𝑐_𝑖𝑥^𝑘_𝑙,0. (19) Together with the definination of Φ_𝑙(𝑥^𝑘_𝑙), implies that (18) holds.
(B2) Since (𝑑^𝑘_𝑙,𝑧^𝑘_𝑙) is the solution of (5), there exists 𝑣∈𝑅 and 𝜇∈𝑅^𝑚𝑙 such that the KKT condition of (5), while 𝑑^𝑘_𝑙=0, we have 𝑣𝑔^𝑘_𝑙+∇𝑐_𝑙𝑥^𝑘_𝑙^𝑇𝜇^𝑘_𝑙=0,1−𝑣−^𝑚𝑙_𝑖=1𝑢^𝑘_𝑖=0,−𝑧^𝑘_𝑙≤0,𝑐_𝑖𝑥^𝑘_𝑙−𝑧^𝑘_𝑙≤0,𝑖∈𝐼_𝑙,𝑣𝑧^𝑘_𝑙𝑢=0,𝑣≥0,_𝑘𝑖𝑐_𝑖𝑥^𝑘_𝑙−𝑧^𝑘_𝑙=0,𝑢^𝑘_𝑖≥0,𝑖∈𝐼_𝑙. (20)(21)(22)(23)(24) By the definition of Φ_𝑙(𝑥^𝑘_𝑙) and (22), Φ_𝑙(𝑥^𝑘_𝑙)≤𝑧^𝑘_𝑙. Then, from (18) and 𝑑^𝑘_𝑙=0, Φ_𝑙(𝑥^𝑘_𝑙)≥𝑧^𝑘_𝑙. From (24), it follows that 𝑢^𝑘_𝑖=0,∀𝑖∉𝐼⁰_𝑙=𝑖∈𝐼_𝑙∶𝑐_𝑖𝑥^𝑘_𝑙=Φ_𝑙𝑥^𝑘_𝑙. (25) Hence, it follows from Assumption (A3) Together with (20) and (21) that 𝑣>0. For (23) we have 𝑧^𝑘_𝑙=0.
From (22) and (24), we have 𝑢^𝑘_𝑖𝑐_𝑖𝑥^𝑘_𝑙=0,𝑢^𝑘_𝑖≥0,𝑖∈𝐼_𝑙,𝑐_𝑖𝑥^𝑘_𝑙≤0,𝑖∈𝐼_𝑙, (26) due to (20) and set 𝜌^𝑘_𝑖=𝑢^𝑘_𝑖/𝑣, we have 𝑔^𝑘_𝑙+∇𝑐_𝑙𝑥^𝑘_𝑙^𝑇𝜌^𝑘_𝑙=0. (27)
Taking into account separability of constraints 𝑐(𝑥), and let 𝐼=∪^𝑝_𝑙=1𝐼_𝑙, then 𝜌^𝑘_𝑖𝑐_𝑖𝑥_𝑘=0,𝜌^𝑘_𝑖𝑐≥0,𝑖∈𝐼,_𝑖𝑥_𝑘𝑔≤0,𝑖∈𝐼,_𝑘+∇𝑐_𝑖𝑥_𝑘^𝑇𝜌_𝑘=0; (28)
that is, 𝑥_𝑘 is a KKT point of (1).

Lemma 6 (see [17, Lemma 3]). In Step 3, the line search procedure is well defined.

Lemma 7 (see [17, Lemma 4]). Under Assumption (A5), the Restoration Algorithm (Algorithm B) terminates finitely.

From Lemmas 6 and 7, we know that the Algorithm A and Algorithm B are well implemented.

The following Lemma implies that every cluster point of {𝑥_𝑘} which generated by Algorithm A is a feasible point of problem (1).

Lemma 8. Let {𝑥_𝑘} be an infinite sequence generated by Algorithm A, then ℎ(𝑥_𝑘)→0(𝑘→∞).

Theorem 9. Suppose {𝑥_𝑘} is an infinite sequence generated by Algorithm A, 𝑑^𝑘_𝑙,𝑙=1,…,𝑝 is the solution of (5), 𝑑_𝑘=(𝑑^𝑘₁,…,𝑑^𝑘_𝑝), then lim_𝑘→∞‖𝑑^𝑘_𝑙‖=0, 𝑙=1,…,𝑝.

Proof. By Assumption (A1), there exists a point 𝑥^∗ such that 𝑥_𝑘→𝑥^∗ for 𝑘∈𝐾, where 𝐾 is an infinite index set. By Algorithm A and Lemma 8, we consider the following two possible cases.

Case I. 𝐾₀={𝑘∈𝐾∣𝑔^𝑇_𝑘𝑑_𝑘≤−(1/2)𝑑^𝑇_𝑘𝐻_𝑘𝑑_𝑘} is an infinite index set. In this case, we have 𝑓𝑥_𝑘+𝜆_𝑘𝑑_𝑘≤𝑓𝑥_𝑙(𝑘)−𝜆_𝑘𝜇𝑑^𝑇_𝑘𝐻_𝑘𝑑_𝑘=𝑓𝑥_𝑙(𝑘)−𝜆_𝑘𝜇^𝑝_𝑙=1𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙. (29)

Since 𝑚(𝑘+1)≤𝑚(𝑘)+1, we obtain 𝑓𝑥_𝑙(𝑘+1)=max_{()0≤𝑗≤𝑚𝑘+1}𝑓𝑥_{𝑘+1−𝑗}≤max_{(𝑘)0≤𝑗≤𝑚+1}𝑓𝑥_{𝑘+1−𝑗}=max𝑓𝑥_𝑙(𝑘),𝑓𝑥_𝑘+1=𝑓𝑥_𝑙(𝑘). (30)

Hence for 𝑚(𝑘)≤𝑀, it holds 𝑓𝑥_𝑙(𝑘)≤𝑓𝑥_{𝑙(𝑙(𝑘)−1)}−𝜆_{𝑙(𝑘)−1}𝜇𝑑^𝑇_{𝑙(𝑘)−1}𝐻_{𝑙(𝑘)−1}𝑑_{𝑙(𝑘)−1}=𝑓𝑥_{𝑙(𝑙(𝑘)−1)}−𝜆_{𝑙(𝑘)−1}𝜇^𝑝_𝑙=1𝑑_𝑙^{𝑙(𝑘)−1𝑇}𝐻_𝑙^{𝑙(𝑘)−1}𝑑_𝑙^{𝑙(𝑘)−1}. (31)

Since 𝑓 is bounded below, {𝑓(𝑥_𝑙(𝑘))} converges. Due to (31), we can obtain lim_𝑘→∞𝜆_{𝑙(𝑘)−1}𝜇^𝑝_𝑙=1𝑑_𝑙^{𝑙(𝑘)−1𝑇}𝐻_𝑙^{𝑙(𝑘)−1}𝑑_𝑙^{𝑙(𝑘)−1}=0. (32)

By Lemma 6, there exists 𝜆>0 such that 𝜆_{𝑙(𝑘)−1}≥𝜆. By Assumption (A4), we obtain lim_𝑘→∞𝑑_𝑙^{𝑙(𝑘)−1}=0,𝑙=1,…,𝑝. (33)

From the uniform continuity of 𝑓(𝑥), this implies that lim_𝑘→∞𝑓𝑥_{𝑙(𝑘)−1}=lim_𝑘→∞𝑓𝑥_𝑙(𝑘). (34)

Let 𝑙(𝑘)=𝑙(𝑘+𝑀+2), for any given 𝑗≥1, we can have lim_𝑘→∞𝑑_𝑙^{𝑙(𝑘)−𝑗}=0,𝑙=1,…,𝑝,lim_𝑘→∞𝑓𝑥_{𝑙(𝑘)−𝑗}=lim_𝑘→∞𝑓𝑥_𝑙(𝑘). (35)(36)

If 𝑗=1, since {𝑙(𝑘)}⊂{𝑙(𝑘)}, (35) and (36) follow from (33) and (34). For a given number 𝑗>1, we have 𝑓𝑥_{𝑙(𝑘)−𝑗}≤𝑓𝑥_{𝑙(𝑙(𝑘)−𝑗−1)}−𝜆_{𝑙(𝑘)−𝑗−1}𝜇𝑑^𝑇_{𝑙(𝑘)−𝑗−1}𝐻_{𝑙(𝑘)−𝑗−1}𝑑_{𝑙(𝑘)−𝑗−1}. (37)

Using the same arguments above, we obtain lim_𝑘→∞𝑑_𝑙^{𝑙(𝑘)−𝑗−1}=0,𝑙=1,…,𝑝,lim_𝑘→∞𝑓𝑥_{𝑙(𝑘)−𝑗−1}=lim_𝑘→∞𝑓𝑥_𝑙(𝑘). (38)

Therefore (35) and (36) hold for any given 𝑗≥1.

Now for any 𝑘, 𝑥_𝑘+1=𝑥_𝑙(𝑘)−∑^{𝑙(𝑘)−𝑘−1}_𝑗=1𝜆_{𝑙(𝑘)−𝑗}𝑑_{𝑙(𝑘)−𝑗}, from Remark 1, we obtain 𝑙(𝑘)=𝑙(𝑘+𝑀+2)<𝑘+𝑀+2, that is 𝑙(𝑘)−𝑘−1≤𝑀+1. Since 𝑑_{𝑙(𝑘)−𝑗}={𝑑^l₁^{(𝑘)−𝑗},…,𝑑_𝑝^{𝑙(𝑘)−𝑗}} together with (35), we can obtain lim_𝑘→∞‖𝑑_{𝑙(𝑘)−𝑗}‖=0, 1≤𝑗≤𝑙(𝑘)−𝑘−1. Therefore, lim_𝑘→∞𝑥_𝑘+1−𝑥_𝑙(𝑘)=0. (39)

Since {𝑓(𝑥_𝑙(𝑘))} admits a limit, by the uniform continuity of 𝑓 on 𝑆, it holds lim_𝑘→∞𝑓𝑥_𝑘=lim_𝑘→∞𝑓𝑥_{𝑙(𝑘)−1}=lim_𝑘→∞𝑓𝑥_𝑙(𝑘). (40)

Then by the relation (29), we have lim^𝑝_{𝑘→∞𝑙=1}𝜆_𝑘𝜇𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙=0. (41)

Using the same arguments for deriving (33), we obtain that lim_𝑘→∞𝑑^𝑘_𝑙=0,𝑙=1,…,𝑝. (42)

Case II. 𝐾₀ is a finite index set, which implies 𝐾₁={𝑘∈𝐾∣𝑔^𝑇_𝑘𝑑_𝑘>−(1/2)𝑑^𝑇_𝑘𝐻_𝑘𝑑_𝑘} is an infinite index set.

If lim_𝑘→∞‖𝑑^𝑘_𝑙‖=0,𝑙=1,…,𝑝 does not hold, then there exists a positive number 𝜀>0 and an infinite index set 𝐾₂, such that ‖𝑑^𝑘_𝑙‖>𝜀, for all 𝑘∈𝐾₂⊂𝐾₁. Since (𝑑^𝑘_𝑙,𝑧^𝑘_𝑙) is the solution of 𝑄𝑃_𝑙. By the KKT condition of (5), we have 𝐻^𝑘_𝑙𝑑^𝑘_𝑙+𝑣𝑔^𝑘_𝑙+∇𝑐_𝑙𝑥^𝑘_𝑙^𝑇𝜇^𝑘_𝑙=0,1−𝑣−^𝑚𝑙_𝑖=1𝑢_𝑖𝑔=0,_𝑙𝑘𝑇𝑑^𝑘_𝑙−𝑧^𝑘_𝑙𝑢≤0,_𝑖𝑐_𝑖𝑥^𝑘_𝑙+∇𝑐_𝑖𝑥^𝑘_𝑙^𝑇𝑑^𝑘_𝑙−𝑧^𝑘_𝑙=0,𝑢_𝑖≥0,𝑖∈I_𝑙. (43)(44)(45)(46)

Since 𝑣≥0, 𝑢_𝑖≥0 and (44), we obtain 0≤𝑢_𝑖≤1, therefore there exists 𝑀_𝑙 satisfied ‖𝑢^𝑘_𝑙‖≤𝑀_𝑙.

By Lemma 8, ℎ(𝑥_𝑘)→0 implies ℎ_𝑙(𝑥^𝑘_𝑙)→0. Hence there exists 𝑘₀, such that for 𝑘>𝑘₀,𝑘∈𝐾₂, it holds ℎ_𝑙𝑥^𝑘_𝑙≤𝑎𝜀²2𝑀_𝑙≤𝑎𝑑^𝑘_𝑙²2𝑀_𝑙≤𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙2𝑀_𝑙,𝑙=1,…,𝑝. (47)

Consequently, from (43), we deduce 𝑣𝑔_𝑙^𝑘𝑇𝑑^𝑘_𝑙=−𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙−𝜇_𝑙^𝑘𝑇∇𝑐_𝑙𝑥^𝑘_𝑙𝑑^𝑘_𝑙=−𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙+^𝑚𝑙_𝑖=1𝑢_𝑖𝑐_𝑖𝑥^𝑘_𝑙−𝑧^𝑘_𝑙≤−𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙+𝑀_𝑙ℎ_𝑙𝑥^𝑘_𝑙−^𝑚𝑙_𝑖=1𝑢_𝑖𝑧^𝑘_𝑙≤−𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙+𝑀_𝑙ℎ_𝑙𝑥^𝑘_𝑙−^𝑚𝑙_𝑖=1𝑢_𝑖𝑔_𝑙^𝑘𝑇𝑑^𝑘_𝑙. (48)

Together with (44), (47), and transpose 𝑔_𝑙^𝑘𝑇𝑑^𝑘_𝑙=^𝑚𝑙_𝑖=1𝑢_𝑖+𝑣𝑔_𝑙^𝑘𝑇𝑑^𝑘_𝑙≤−𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙+𝑀_𝑙ℎ_𝑙𝑥^𝑘_𝑙≤−𝑑^𝑘𝑇_l𝐻^𝑘_𝑙𝑑^𝑘_𝑙+12𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙1=−2𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙. (49)

Taking into account separability of constraints 𝑐(𝑥), then 𝑔^𝑇_𝑘𝑑_𝑘=^𝑝_𝑙=1𝑔_𝑙^𝑘𝑇𝑑^𝑘_𝑙≤−^𝑝_𝑙=112𝑑_𝑙^𝑘𝑇𝐻^𝑘_𝑙𝑑^𝑘_𝑙1=−2𝑑^𝑇_𝑘𝐻_𝑘𝑑_𝑘, (50) which contradicts the definition of 𝐾₁. The proof is complete.

Theorem 10. Suppose {𝑥_𝑘} is an infinite sequence generated by Algorithm A and the assumptions in Theorem 9 hold. Then every cluster point of {𝑥_𝑘} is a KKT point of problem (1).

Proof. By Assumption (A2), there exists 𝑥^∗, such that 𝑥_𝑘→𝑥^∗,𝑘∈𝐾. From Lemma 8, we have ℎ(𝑥_𝑘)→0,𝑘∈𝐾, which means that 𝑥^∗ is a feasible point. From Theorem 9, we have lim_𝑘→∞‖𝑑^𝑘_𝑙‖=0,𝑙=1,…,𝑝. By Lemma 5, 𝑑^𝑘_𝑙→0 implies that 𝑧^𝑘_𝑙→0,𝑙=1,…,𝑝, that is, 𝑑^∗_𝑙=0,𝑙=1,…,𝑝,𝑧^∗_𝑙=0,𝑙=1,…,𝑝.
Compare the KKT condition of (5) with the KKT condition of problem (1), and let 𝜌^∗=(𝜌^∗₁,…𝜌^∗_𝑝),𝜌^∗_𝑙=(𝑢^∗_𝑖/𝑣,𝑖∈𝐼_𝑙),𝐼=∪^𝑝_𝑙=1𝐼_𝑙.
Taking into account separability of constraints 𝑐(𝑥) and 𝑥, we conclude 𝑔^∗+𝜌^∗∇𝑐^∗=0,𝑐𝑥^∗≤0,𝜌^∗_𝑖𝜌≥0,_𝑖∗𝑇𝑐_𝑖𝑥^∗=0,𝑖∈𝐼. (51)
Therefore 𝑥^∗ is a KKT point of problem (1).

4. Numerical Results

To get some insight into computational properties of our approach in Section 2, we considered the same test problems taken from [1]. Choose Problem 1 of [1] as follows: min^𝑝−1𝑝_𝑖=1^𝑣_{𝑗=𝑖+1𝑡=1}𝑥_{(𝑖−1)𝑣+𝑡}𝑥_{(𝑗−1)𝑣+𝑡}s.t.^𝑣_𝑡=1𝑥²_{()𝑖−1𝑣+𝑡}−1=0,𝑖=1,…,𝑝. (52)

In our test, we only verify our convergence theory without comparing with serial algorithms. In the future, we will propose the convergent rate of the parallel algorithm.

Not having access to a parallel computer, we have carried out on Graphics Processing Unit (GPU) to solve them and implement the algorithm on Compute Unified Device Architecture (CUDA) [20]. However, the input and output are implemented by CPU. And all the subproblems are solved on blocks which are constructed into many threads. If there are 𝑝 blocks in GPU, which is equivalent to 𝑝 processors in a parallel computer. Many threads of the block can achieve a large number of vector calculations in current block.

We have implemented Algorithm A in C language. All codes were run under Visual Studio 2008 and Cuda 4.0 on the DELL Precision T7400. In Table 1, we report the number of iterations and the running times for Algorithms A on Problem (52). The results in Table 1 confirm that the proposed approach certainly makes sense. We solve the subproblem (5) which is a positive semidefinite quadratic programs by the modified interior-point methods. And the Restoration Algorithm of [17] is solved by Trust-region Approach combined with local pattern search methods. All algorithms are coded by C language without using any function from optimization Toolbox; hence the results are not the best. In Table 1, 𝑛 denotes the dimensions of test problem and 𝜈 denotes the variable number of one block. There is less iteration number of 𝑄𝑃_𝑙 subproblem with more numbers of blocks for high-dimensional problem, which shows the algorithm is reliable.

Table 1 

Results for Algorithm A.

5. Conclusion

The PVD algorithm which was proposed in 1994 is used to solve unconstrained optimization problems or has a special case of convex block-separable constraints. In 1997, M. V. Solodov proposed useful generalizations that consist, for the general unconstrained case, of replacing exact global solution of the subproblems by a certain natural sufficient descent condition, and, for the convex case, of inexact subproblem solution in the PVD algorithm. M. V. Solodov proposed a PVD approach in 1998 which applied to problems with general convex constraints directly and show that the algorithm converges. In 2002, C. A. Sagastizábal et al. propose two variants of PVD for the constrained case: without assuming convexity of constraints, but assuming block-separable structure and for inseparable constraints, but assuming convexity. In this paper, we propose the modified algorithm of [1] used to the general constraints but with block-separable structure.

In a word, the algorithm above is a special structure of the objective function or constraints with a special structure. In the further, we will study the parallel algorithm with general inseparable constraints.

Acknowledgments

This paper was supported by National Natural Science Foundation of China (11101420, 71271204) and other foundations (2010BSE06047). The authors would like to thank Dr. Jiang Zhipeng for valuable work on numerical experiments.