# Sufficient Descent Polak-Ribière-Polyak Conjugate Gradient Algorithm for Large-Scale Box-Constrained Optimization

Authors: Qiuyu Wang; Yingtao Che
Journal: Abstract and Applied Analysis
 (2014)

---

## Abstract

A practical algorithm for solving large-scale box-constrained optimization problems is developed, analyzed, and tested. In the proposed algorithm, an identification strategy is involved to estimate the active set at per-iteration. The components of inactive variables are determined by the steepest
descent method at first finite number of steps and then by conjugate gradient method subsequently. Under some appropriate conditions, we show that the algorithm converges globally. Numerical experiments and comparisons by using some box-constrained problems from CUTEr library are reported. Numerical comparisons illustrate that the proposed method is promising and competitive with the well-known method—L-BFGS-B.


---

## Body

## 1. Introduction

We consider the box-constrained optimization problem(1)
min
⁡
x
∈
ℝ
n
f
(
x
)
,
l
≤
x
≤
u
,
where f
:
ℝ
n
→
ℝ is continuously differentiable, and l, u
∈
ℝ
n with l
<
u. The gradient of f at x
k is ∇
f
(
x
k
)
∈
ℝ
n, and its ith component is ∇
f
i
(
x
k
) (∈
ℝ) or g
i
k for the sake of simplicity. We define the feasible region of (1) as Ω; that is,
(2)
Ω
=
{
x
∈
ℝ
n
:
l
≤
x
≤
u
}
.
We say that a vector x
¯
∈
Ω is a stationary point for problem (1) if it satisfies
(3)
l
i
=
x
¯
i
⟹
∇
f
i
(
x
¯
)
≥
0
,
l
i
<
x
¯
i
<
u
i
⟹
∇
f
i
(
x
¯
)
=
0
,
x
¯
i
=
u
i
⟹
∇
f
i
(
x
¯
)
≤
0
.
Problem (1) is very important in practical optimization, because numerous practical problems can be converted into this form. In addition, problem (1) is often treated as a subproblem of augmented Lagrangian or penalty schemes for general constrained optimization. Hence, the development of numerical algorithms to efficiently solve (1), especially when the dimension is large, is important both in theory and application. Moreover, the box-constrained optimization problems have been received much attention in recent decades. We refer to the excellent paper [1] for a good review.Many algorithms for solving this type of problems are based on active set strategies, for example, [2]. In this class of methods, a working set is used to estimate the set of active constraints at the solution and it is updated at per-iteration. The early active set methods are quite efficient for relatively lower dimensional problems but are typically unattractive for large-scale ones [3]. The main reason is that at most one constraint can be added to or dropped from the active set at each step. However, this slows down the rate of the convergence. Based on this consideration, more and more scholars are engaged in designing active set methods aiming to make rapid changes against/from incorrect predictions, for example, [3–8].LetΩ
=
ℝ
n. It is clear that problem (1) reduces to the following unconstrained optimization problem:
(4)
min
⁡
x
∈
ℝ
n
f
(
x
)
.
It can be solved by almost any existing effective methods such as the nonlinear conjugate gradient method, limited memory BFGS method, and spectral gradient method. These methods are much suitable for solving (4) with n being assumed to be large due to their simplicity and low storage requirements.In conjugate gradient methods scheme, Polak-Ribière-Polyak (PRP) method is generally believed to be the most efficient from computation point of view. But its global convergence for nonconvex nonlinear function is still uncertain [9]. To overcome this shortcoming, various techniques have been proposed. In particular, Zhang et al. [10] introduced a modified MPRP method in which the generated direction is descent independent of line search rule [11]. Global convergence is guaranteed by means of a variation of the Armijo-like line search strategy of Grippo and Lucidi [12]. Moreover, MPRP is also considered to be competent to solve the box-constrained optimization problem (1) based on projected gradient techniques [13]. However, to the best of our knowledge, papers on nonlinear conjugate gradient for solving box-constrained minimization problems are relatively fewer.In this paper, from a different point of view, we further study the applications of the MPRP method in [10] for solving box-constrained minimization problems. First, motivated by the work of Facchinei et al. [3], we estimate a lower dimensional free subspace at per-iteration, which has the potencies of being numerically cheaper when applied to solve large-scale problems. The search direction at per-iteration consists of two parts: some of the components are defined simply; the others are determined by the nonlinear gradient method. The remarkable feather of the proposed method is that all the generated points are feasible without requiring gradient projection. Under some appropriate conditions, we show that the method converges globally. We present experimental results comparing the developed method with the well-known algorithm L-BFGS-B.The paper is organized as follows. For easy comprehension of the proposed algorithm, we briefly recall an active set identification technique and then state the steps of our new algorithm in Section2. In Section 3, we show that the proposed method converges globally. In Section 4, we test the performance of the algorithm and do some numerical comparisons. Finally, we draw some conclusions in Section 5. In the rest of this paper, the symbol ∥
·
∥ denotes the Euclidean norm of vectors.


## 2. Motivation and Algorithm

In this section, we introduce a new algorithm which is defined by(5)
x
k
+
1
=
x
k
+
α
k
d
k
,
where d
k is a search direction and α
k is the corresponding step length. To construct our algorithm, throughout this paper, we assume that the following assumptions hold.Assumption 1.
The level setℱ
=
{
x
∈
ℝ
n
:
f
(
x
)
≤
f
(
x
0
)
}
∩
Ω is compact.Assumption 2.
In some neighborhood𝒩 of ℱ, the gradient ∇
f
(
x
) is Lipschitz continuous; that is, there exists a constant ℒ
>
0 such that
(6)
∥
∇
f
(
x
)
-
∇
f
(
y
)
∥
≤
ℒ
∥
x
-
y
∥
,
∀
x
,
y
∈
𝒩
.Assumption 3.
Strict complementarity condition holds atx
¯; that is, the strict inequalities hold in the first and the last implications of (3).In [3], Facchinei et al. presented an active set Newton algorithm for solving problem (1). Their algorithm possesses some favorable properties, such as fast local convergence and feasibility of all iterations. In addition, only a lower dimensional quadratic program subproblem needs to be solved at per-iteration. We now simply recall the active set identification technique in [3]. Let a
i
(
x
) and b
i
(
x
) be nonnegative continuous and bounded functions defined on Ω, such that if x
i
=
l
i or x
i
=
u
i, then a
i
(
x
)
>
0 or b
i
(
x
)
>
0, respectively. For any x
∈
Ω, define the index set L
(
x
), F
(
x
) and U
(
x
) as
(7)
L
(
x
)
=
{
i
:
x
i
≤
l
i
+
a
i
(
x
)
∇
f
i
(
x
)
}
,
U
(
x
)
=
{
i
:
x
i
≥
u
i
+
b
i
(
x
)
∇
f
i
(
x
)
}
,
F
(
x
)
=
{
1
,
…
,
n
}
∖
(
L
∪
U
)
.
The set L
(
x
)
∪
U
(
x
) is an estimate of the active set at point x. Facchinei et al. [3] showed that with Assumption 3, when x is sufficiently close to x
¯, the index set L
(
x
)
∪
U
(
x
) is an exact identification of L
(
x
¯
)
∪
U
(
x
¯
). This property will enable us to find the active set at the stationary point x
¯ in a finite number of steps. This active set identification technique (7) has motivated much additional studies on box-constrained optimization problems, for example, [14–18].Using the active set estimation technique in (7), we deduce the search direction used in our algorithm in detail now. For the sake of simplicity, let L
k
=
L
(
x
k
),  U
k
=
U
(
x
k
),  F
k
=
F
(
x
k
), and let |
·
| be the number of elements in the corresponding vector let Z
k be the matrix whose columns are {
e
i
:
i
∈
F
k
}, where e
i is the ith column of the identity matrix. Additionally, we denote, respectively, the ith component of x
k and d
k by x
i
k and d
i
k.Letd
k
∈
ℝ
n, in which
(8)
d
i
k
=
l
i
-
x
i
k
,
i
∈
L
k
,
d
i
k
=
u
i
-
x
i
k
,
i
∈
U
k
.
Now, we restrict our attention to define components with index i
∈
F
k. Let
(9)
d
~
k
=
{
-
g
~
0
,
if
k
=
0
,
-
g
~
k
,
if
k
>
0
,
F
k
≠
F
k
-
1
,
-
g
~
k
+
β
k
d
~
k
-
1
-
θ
k
y
~
k
-
1
,
if
k
>
0
,
F
k
=
F
k
-
1
,
where
(10)
β
k
=
(
g
~
k
)
⊤
y
~
k
-
1
∥
g
~
k
-
1
∥
2
,
θ
k
=
(
g
~
k
)
⊤
d
~
k
-
1
∥
g
~
k
-
1
∥
2
,
in which g
~
k
=
Z
k
⊤
g
k, g
~
k
-
1
=
Z
k
⊤
g
k
-
1, and y
~
k
-
1
=
g
~
k
-
g
~
k
-
1. Clearly,  d
~
k
∈
ℝ
|
F
k
| and therefore Z
k
d
~
k
∈
ℝ
n. It is not difficult to see that the definition of the search direction d
~
k ensures the sufficient descent property; that is
(11)
(
g
~
k
)
⊤
d
~
k
=
-
∥
g
~
k
∥
2
.
The definition of the search direction (9) was designed by Zhang et al. [10] in full space ℝ
n. Here, it is stated for slightly different circumstances, but it is easy to verify that it is still valid in our context.In order to guarantee the feasibility of all iterations and the descent of the objective function at each iteration, we define the positive scalarξ
k
* by
(12)
ξ
k
*
=
max
⁡
{
ξ
∣
ξ
≤
1
,
l
i
-
x
i
k
≤
ξ
(
Z
k
d
~
k
)
i
≤
u
i
-
x
i
k
,
i
∈
F
k
}
.
Note that when the strict complementarity condition holds and l
<
x
<
u, then ξ
k
* is always well defined unless the stationary point is achieved. Combining with (8) and (9), the compact form of d
k can be rewritten as
(13)
d
i
k
=
{
l
i
-
x
i
k
,
i
∈
L
k
,
u
i
-
x
i
k
,
i
∈
U
k
,
ξ
k
*
(
Z
k
d
~
k
)
i
,
i
∈
F
k
.Remark 4.
Observe that the search direction of inactive variables in (9) contain both gradient and direction information at the current or the previous steps, while the set of inactive variables might be changed as the iteration progresses. In particular, the sets F
k
-
1 and F
k might not be the same. In this situation, we set d
~
k
-
1
=
-
g
~
k
-
1, that is, a negative gradient direction. Otherwise, we use the conjugate gradient formula to update the search direction in the free subspace.Remark 5.
Considering the case when the current iterationx
k is still not an optimal point, the value of g
~
k might be sufficiently small. In this case, for practical purpose, it might be efficient to restrict ∥
g
~
k
∥ bounded from zero.To list our algorithm in detail, we first show some useful properties of the direction. The following result shows that wheneverd
k
≠
0, it is a descent direction for objective function f at current point x
k. The property is very important to our algorithm.Lemma 6.
Ifd
k is attained by (13), it satisfies
(14)
(
g
k
)
⊤
d
k
≤
0
,
and the equality holds if and only if d
k
=
0.Proof.
By the definition of search directiond
k, we have
(15)
(
d
k
)
⊤
g
k
=
∑
i
∈
L
k
‍
(
l
i
-
x
i
k
)
g
i
k
+
∑
i
∈
U
k
‍
(
u
i
-
x
i
k
)
g
i
k
-
∥
Z
k
⊤
g
k
∥
2
≤
0
.
The above relation can be easily obtained from the definition of the direction d
k in (13). Notice that the strict complementary condition holds, which shows (
d
k
)
⊤
g
k
=
0 if and only if d
k
=
0.As an immediate consequence of the definition of the search direction, we have the following fact.Lemma 7.
Letd
k be determined by (13) and assume that d
k
≠
0; then,
(16)
x
k
+
d
k
∈
Ω
.The above lemma shows that for all0
≤
α
≤
1 such that
(17)
P
Ω
(
x
k
+
α
d
k
)
=
x
k
+
α
d
k
.
Since the search direction d
k has been computed, a step length α
k
∈
(
0,1
] ought to be determined in order to obtain the next iteration as (5). In this paper, we consider a backtracking line search procedure, where a decreasing sequence of {
ρ
j
} is tried until the smallest j
k is found which satisfies
(18)
f
(
x
k
+
ρ
j
d
k
)
≤
f
(
x
k
)
-
δ
∥
ρ
j
d
k
∥
2
,
where ρ
∈
(
0,1
), δ
∈
(
0,1
/
2
). Then a step length is accepted corresponding to α
k
=
ρ
j
k.We are now ready to formally state the overall algorithm for solving the box-constrained optimization problems (1) as follows.Algorithm 8 (SDPRP).
The steps of the algorithm is given as follows.
Step  0.  Given starting point x
0
∈
Ω, constants δ
∈
(
0,1
/
2
), ρ
∈
(
0,1
) and g
max
⁡
>
g
min
⁡
>
0. Set k
:
=
0.
Step  1.  Determine L
k, U
k, and F
k according to (7).
Step  2.  If F
k
≠
F
k
-
1, set d
~
k
=
-
g
~
k. Otherwise, set
(19)
∥
g
~
k
-
1
∥
2
=
min
⁡
{
g
max
⁡
,
max
⁡
{
g
min
⁡
,
∥
g
~
k
-
1
∥
2
}
}
.
Determine d
k by (13).
Step  3.  If d
k
=
0, then stop.
Step  4.  Find the smallest integer j (say j
k) satisfying (18). Let α
k
=
ρ
j
k.
Step  5.  Set x
k
+
1
:
=
x
k
+
α
k
d
k.
Step  6.  Let k
:
=
k
+
1, go to Step  1.Remark 9.
Now suppose that the active sets have been identified after a finite number of steps; then, we haved
i
k
=
0 for all i
∈
L
k
∪
U
k and d
i
k
=
ξ
k
*
(
Z
k
d
~
k
)
i for all i
∈
F
k. If α
k
≠
1, by the line search process, we know that α
k
′
=
ρ
-
1
α
k does not satisfy (18). That is,
(20)
f
(
x
k
+
ρ
-
1
α
k
d
k
)
-
f
(
x
k
)
>
-
δ
ρ
-
2
α
k
2
∥
d
k
∥
2
.
By the mean-value theorem, there is aζ
k
∈
(
0,1
) such that
(21)
f
(
x
k
+
ρ
-
1
α
k
d
k
)
-
f
(
x
k
)
=
ρ
-
1
α
k
g
(
x
k
+
ζ
k
ρ
-
1
α
k
d
k
)
⊤
d
k
=
ρ
-
1
α
k
(
g
k
)
⊤
d
k
+
ρ
-
1
α
k
[
g
(
x
k
+
ζ
k
ρ
-
1
α
k
d
k
)
-
g
k
]
⊤
d
k
≤
ρ
-
1
α
k
(
g
k
)
⊤
d
k
+
ℒ
ρ
-
2
α
k
2
∥
d
k
∥
2
,
where ℒ is the Lipschitz constant. It follows from [10, Lemma 3.1] that there exists a positive constant M such that ∥
d
k
∥
≤
M. Substituting the above inequality into (20), we have
(22)
α
k
≥
ρ
ℒ
+
δ
(
g
k
)
⊤
d
k
∥
d
k
∥
2
=
ρ
ℒ
+
δ
∥
g
~
k
∥
2
∥
d
~
k
∥
2
≥
ρ
M
(
ℒ
+
δ
)
∥
g
~
k
∥
2
.
It follows that when ∥
g
~
k
∥
2
≠
0, the line search step is well defined, so is the whole algorithm.


## 3. Convergence Analysis

In this section, we show that Algorithm8. converges globally. The following lemma gives a sufficient and necessary condition for the global convergence of Algorithm 8.Lemma 10.
Let{
x
k
} be a sequence of iterations generated by Algorithm 8, and let d
k be a search direction defined by (13). Then x
k is a stationary point of (1) if and only if d
k
=
0.Proof.
Letd
k
=
0. If i
∈
L
k, then we have
(23)
0
=
d
i
k
=
l
i
-
x
i
k
≥
-
a
i
(
x
k
)
∇
f
i
(
x
k
)
.
Since x
i
k
=
l
i and a
i
(
x
k
)
>
0, the right inequality implies ∇
f
i
(
x
k
)
≥
0. If ∈
U
k, we will prove it similarly. And if ∈
F
k, we can get from (11) that ∇
f
i
(
x
k
)
=
0. Now suppose that x
k is a stationary point of (1); it gets from (3) and (7) that
(24)
L
k
=
{
i
:
x
i
k
=
l
i
}
,
F
k
=
{
i
:
l
i
<
x
i
k
<
u
i
}
,
U
k
=
{
i
:
x
i
k
=
u
i
}
.
Then it follows from (9) and (11) that d
L
k
k
=
0, d
U
k
k
=
0 and d
F
k
k
=
0. Therefore d
k
=
0.Lemma 11.
Let{
x
k
} be a sequence of iterations generated by Algorithm 8 and let d
k be a search direction defined by (13). Suppose that there are subsequences {
x
k
}
𝒦
→
x
¯ and {
d
k
}
𝒦
→
0 as k
→
∞. Then x
¯ is a stationary point of problem (1).Proof.
Since the number of distinct setsL
k, U
k, and F
k is finite and Assumptions 1–3 hold, Assumption 3 ensures that there exists k
0 such that for all k
>
k
0 [3],
(25)
L
k
=
L
¯
,
U
k
=
U
¯
,
F
k
=
F
¯
.
Since {
x
k
}
⊆
ℱ, we obviously have
(26)
l
≤
x
¯
≤
u
.
Furthermore, the fact that {
d
L
k
}
𝒦
→
0 and {
d
U
k
}
𝒦
→
0 implies
(27)
x
¯
L
=
l
L
,
∇
f
L
(
x
¯
)
≥
0
,
x
¯
U
=
u
U
,
∇
f
U
(
x
¯
)
≤
0
.
The proof of ∇
f
F
¯
(
x
¯
)
=
0 can be obtained once we notice that (11) and d
F
k
k
→
0 as k
→
∞.The similar proof of Lemmas10 and 11 can also be found in [3]. To end of this section, now we are ready to establish global convergence of Algorithm 8.Theorem 12.
Suppose that Assumptions1–3 hold. Then the sequence {
x
k
} generated by Algorithm 8 has at least a limit point, and every limit point of this sequence is a stationary point of problem (1).Proof.
It easily follows from Lemma7 that the sequence {
x
k
} of points generated by the Algorithm 8 is contained in the compact set ℱ. Hence, by Assumption 1, there exists at least a limit point of this sequence.
If the sequence{
x
k
} is finite with last point x
¯, then, by Lemma 10, x
¯ is a stationary point of problem (1). So we assume that the sequence is infinite.
Regarding the first part of proof in Lemma11, L
k, U
k, and F
k are constants starting from some iteration index k
0; that is,
(28)
L
k
=
L
¯
,
U
k
=
U
¯
,
F
k
=
F
¯
.
It is easy to see that(29)
d
i
k
=
l
i
k
-
x
i
k
=
0
,
∀
i
∈
L
k
,
d
i
k
=
u
i
k
-
x
i
k
=
0
,
∀
i
∈
U
k
,
as k
→
∞. Now, it remains to prove d
i
k
→
0 (i
∈
F
k) as k
→
∞. According to the global convergence theorem in [10], we have
(30)
liminf
⁡
k
→
∞
∥
g
~
i
k
∥
=
0
,
∀
i
∈
F
k
.
By the definition of d
~
k in (9), it is equivalent to
(31)
liminf
⁡
k
→
∞
∥
d
~
i
k
∥
=
0
,
∀
i
∈
F
k
.
It follows from (29) and (31) that {
d
k
}
k
∈
𝒦
→
0, which shows our claims by Lemma 11. The proof is complete.


## 4. Numerical Experiments

Now, let us report some numerical results attained by our Sufficient Descent Polak-Ribière-Polyak Algorithm—SDPRP. The algorithm is implemented by Fortran77 code in double precision arithmetic. All experiments are run on a PC with CPU Intel Pentium Dual E2140 1.6 GHz, 512 M bytes of SDRAM memory, and Red Hat Linux 9.03 operating system. Our experiments are performed on a set of the nonlinear box-constrained problems from the CUTEr [19] library that have second derivatives available. Since we are interested in large problems, we refine this selection by considering only problems where the number of variables is at least 50. Altogether, we solve 54 problems. The type of objective function and the character of the problems being tested are listed in Table 1.Table 1
Problem set according to the CUTEr classification.
Number of sets
Objective type
Problem character
Classification
1
Others
Academic
OBR2-AN-*-*
2
Others
Modeling
OBR2-MN-*-*
3
Others
Real application
OBR2-RN-*-*
4
Sum of squares
Academic
SBR2-AN-*-*
5
Sum of squares
Modeling
SBR2-AN-*-*
6
Quadratic
Academic
QBR2-AN-*-*
7
Quadratic
Modeling
QBR2-MN-*-*
8
Quadratic
Real application
QBR2-RN-*-*
9
Others
Modeling
OBI2-MY-*-*
10
Quadratic
Academic
QBR2-AY-*-*In the experiments, for easily comparing with other codes, we use the projected gradient errors to measure the quality of the solutions instead of∥
d
k
∥; that is, we force the iteration stopped when
(32)
∥
P
Ω
(
x
k
-
∇
f
(
x
k
)
)
-
x
k
∥
∞
≤
1
0
-
5
,
where P
Ω
(
·
) is the projected gradient of objective function and ∥
·
∥
∞ denotes the maximum absolute component of a vector. Clearly, the stopping criteria (32) are equivalent to ∥
∇
f
(
x
k
)
∥
∞
≤
1
0
-
5 in our method. We also stop the execution of SDPRP when 10000 iterations or 20000 function evaluations are completed without achieving convergence. We choose the initial parameters g
min
⁡
=
1
0
-
7, g
max
⁡
=
1
0
20, ρ
=
0.29, and δ
=
1
0
-
1. Moreover, we also test our proposed methods with different parameters a
i
(
x
) and b
i
(
x
) to see that a
i
(
x
)
=
b
i
(
x
)
=
1
0
-
6
∥
P
Ω
(
∇
f
(
x
0
)
∥ is the best choice. In order to assess the performance of the SDPRP, we then test the well-known method L-BFGS-B (available at http://www.ece.northwestern.edu/~nocedal/lbfgsb.html) [20, 21] for comparison.When running the codes L-BFGS-B, default values are used for all parameters. For L-BFGS-B, besides (32), the algorithm has another built-in stopping test based on the parameterfactr. It is designed to terminate the run when the change in the objective function f is sufficiently small. The default valuefactr=
1.0
d
+
7 for moderate accuracy. Additionally, it can occur sometimes that the line search cannot make any progress, or the value of Nfg reached a prefixed number (=10000); in these cases, the run is also terminated.The numerical results of the algorithms SDPRP and L-BFGS-B are listed in Tables2 and 3, respectively. We used horizontal lines in both tables to divide the selected problems into 10 classes according to Table 1. The columns in Tables 2 and 3 have the following meanings:Problem: name of the problem;
Dim: dimension of the problem;
Iter: number of iterations;
Nf: number of function evaluations;
Ng: number of gradient evaluations;
Time: CPU time in seconds;
Fv: final function value;
Pgnorm1: maximum-norm of the final gradient projection.Table 2
Performance of SDPRP.
Problem
Dim
Iter
Nf
Time
Fv
Pgnorm1
SINEALI
1000
19
64
0.030
-
0.99860
E
+
05
0.54718
E
-
05
BDEXP
1000
10001
10003
6.030
0.92538
E
-
03
0.13189
E
-
04
EXPLIN
120
76
260
0.010
-
0.72376
E
+
06
0.95740
E
-
05
EXPLIN2
120
59
183
0.000
-
0.72446
E
+
06
0.50532
E
-
05
EXPQUAD
120
4007
20001
0.230
-
0.36260
E
+
07
0.28861
E
-
03
MCCORMCK
5000
32
72
0.130
-
0.45666
E
+
04
0.69662
E
-
05
PROBPENL
500
18
173
0.010
0.39919
E
-
06
0.71380
E
-
05
QRTQUAD
120
949
4240
0.030
-
0.36246
E
+
07
0.89349
E
-
05
S368
100
54
136
0.240
-
0.14025
E
+
03
0.42537
E
-
05
HADAMALS
1024
113
479
0.390
0.30658
E
+
05
0.31465
E
-
05
SCOND1LS
5002
5061
20001
81.330
0.11463
E
+
03
0.10058
E
+
00
CHEBYQAD
50
735
4514
6.440
0.53863
E
-
02
0.98736
E
-
05
HS110
200
1
2
0.000
-
0.99601
E
+
40
0.00000
E
+
00
LINVERSE
1999
1527
1565
2.220
0.68200
E
+
03
0.89735
E
-
05
NONSCOMP
5000
67
128
0.110
0.23525
E
-
10
0.42102
E
-
05
QR3DLS
610
2760
20000
5.530
0.17615
E
+
00
0.29073
E
-
01
DECONVB
61
3353
11830
0.260
0.35269
E
-
07
0.95952
E
-
05
QUDLIN
5000
1
2
0.00
-
0.12500
E
+
10
0.00000
E
+
00
BIGGSB1
5000
10001
11535
10.93
0.21870
E
-
01
0.74163
E
-
04
BQPGABIM
50
45
213
0.00
-
0.37903
E
-
04
0.75229
E
-
05
BQPGASIM
50
57
223
0.00
-
0.55198
E
-
04
0.62637
E
-
05
BQPGAUSS
2003
3945
20001
6.78
-
0.18077
E
+
00
0.98485
E
-
01
CHENHARK
5000
6335
20001
9.56
-
0.19998
E
+
01
0.36782
E
-
04
CVXBQP1
10000
14
15
0.05
0.22502
E
+
07
0.00000
E
+
00
HARKERP2
100
10
47
0.00
-
0.50000
E
+
00
0.20897
E
-
07
JNLBRNG1
10000
4001
7427
38.40
-
0.18057
E
+
00
0.96930
E
-
05
JNLBRNG2
10000
1911
6880
23.35
-
0.41486
E
+
01
0.85547
E
-
05
JNLBRNGA
10000
2182
5955
20.17
-
0.27109
E
+
00
0.94192
E
-
05
JNLBRNGB
10000
4981
20001
54.23
-
0.62068
E
+
01
0.23366
E
-
01
NCVXBQP1
10000
1
2
0.01
-
0.19855
E
+
11
0.27756
E
-
16
NCVXBQP2
10000
3180
20000
16.18
-
0.13340
E
+
11
0.13402
E
+
01
NCVXBQP3
10000
3184
20000
17.11
-
0.65583
E
+
10
0.92309
E
+
01
NOBNDTOR
5476
893
1462
3.54
-
0.44993
E
+
00
0.51489
E
-
05
OBSTCLAE
10000
5261
5794
40.60
0.18865
E
+
01
0.97210
E
-
05
OBSTCLAL
10000
813
1121
6.29
0.18865
E
+
01
0.83060
E
-
05
OBSTCLBL
10000
3509
3872
26.90
0.72722
E
+
01
0.56157
E
-
05
OBSTCLBM
10000
3362
3730
25.97
0.72722
E
+
01
0.87301
E
-
05
OBSTCLBU
10000
2126
2454
16.71
0.72722
E
+
01
0.86796
E
-
05
PENTDI
1000
11
22
0.00
-
0.75000
E
+
00
0.31759
E
-
05
TORSION1
10000
1550
2401
12.51
-
0.42726
E
+
00
0.72463
E
-
05
TORSION2
10000
3227
4745
26.11
-
0.42725
E
+
00
0.90789
E
-
05
TORSION3
10000
394
712
3.23
-
0.12138
E
+
01
0.62159
E
-
05
TORSION4
10000
2183
2957
17.40
-
0.12138
E
+
01
0.76428
E
-
05
TORSION5
10000
114
228
1.26
-
0.28604
E
+
01
0.86795
E
-
05
TORSION6
10000
1639
1847
12.37
-
0.28604
E
+
01
0.54347
E
-
05
TORSIONA
10000
1410
2273
13.26
-
0.41838
E
+
00
0.97648
E
-
05
TORSIONB
10000
2662
4351
25.30
-
0.41838
E
+
00
0.97029
E
-
05
TORSIONC
10000
433
747
4.41
-
0.12044
E
+
01
0.65687
E
-
05
TORSIOND
10000
2269
3019
20.06
-
0.12044
E
+
01
0.69877
E
-
05
TORSIONE
10000
97
194
0.95
-
0.28507
E
+
01
0.68035
E
-
05
TORSIONF
10000
1589
1775
13.42
-
0.28507
E
+
01
0.94166
E
-
05
ODNAMUR
11130
10001
16448
44.78
0.11336
E
+
05
0.21497
E
+
02
GRIDGENA
6218
620
20000
16.58
-
0.19243
E
+
16
0.99274
E
+
01
NOBNDTOR
5476
893
1462
3.69
-
0.44993
E
+
00
0.51489
E
-
05Table 3
Performance of L-BFGS-B.
Problem
Dim
Nf
Ng
Time
Fv
Pgnorm1
SINEALI
1000
44
44
0.03
-
0.9990
E
+
05
0.8973
E
-
03
BDEXP
1000
18
18
0.01
0.4808
E
-
03
0.6140
E
-
05
EXPLIN
120
40
40
0.00
-
0.5150
E
+
08
0.3206
E
+
07
EXPLIN2
120
43
43
0.00
-
0.7833
E
+
08
0.5921
E
+
07
EXPQUAD
120
8
8
0.00
-
0.3163
E
+
07
0.1181
E
+
05
MCCORMCK
5000
10150
10150
40.29
-
0.1221
E
+
12
0.1426
E
+
03
PROBPENL
500
4
4
0.00
0.3992
E
-
06
0.3071
E
-
06
QRTQUAD
120
13349
13349
0.63
-
0.3561
E
+
08
0.3110
E
+
06
S368
100
21
21
0.07
-
0.4084
E
+
02
0.2220
E
+
02
HADAMALS
1024
20
20
0.04
0.3165
E
+
05
0.4539
E
-
05
SCOND1LS
5002
25
25
0.15
0.6374
E
+
05
0.7813
E
+
03
CHEBYQAD
50
21
21
0.15
0.1395
E
-
01
0.1642
E
+
01
HS110
200
21
21
0.02
-
0.1478
E
+
39
0.3285
E
+
37
LINVERSE
1999
195
195
0.32
0.6810
E
+
03
0.4605
E
-
04
NONSCOMP
5000
42
42
0.09
0.7536
E
-
11
0.6216
E
-
05
QR3DLS
610
11054
11054
10.24
0.4538
E
-
02
0.1914
E
-
01
DECONVB
61
291
291
0.02
0.5845
E
-
07
0.8405
E
-
05
QUDLIN
5000
314
314
0.16
-
0.1315
+
301
0.1380
+
151
BIGGSB1
5000
4452
4452
8.36
0.1481
E
-
03
0.8917
E
-
05
BQPGABIM
50
41
41
0.00
-
0.1839
E
-
03
0.7289
E
-
05
BQPGASIM
50
41
41
0.00
-
0.1839
E
-
03
0.7289
E
-
05
BQPGAUSS
2003
9873
9873
11.62
-
0.7222
+
301
0.1179
+
154
CHENHARK
5000
10323
10323
21.54
-
0.8985
E
+
10
0.1722
E
+
03
CVXBQP1
10000
6798
6798
30.27
0.1339
E
-
07
0.5419
E
-
03
HARKERP2
100
11329
11329
1.26
-
0.3318
E
+
02
0.3309
E
+
00
JNLBRNG1
10000
333
333
2.99
-
0.1146
E
+
01
0.9604
E
-
05
JNLBRNG2
10000
998
998
9.02
-
0.1065
E
+
03
0.9553
E
-
05
JNLBRNGA
10000
440
440
3.60
-
0.2179
E
+
01
0.7975
E
-
05
JNLBRNGB
10000
6375
6375
52.41
-
0.9552
E
+
03
0.9473
E
-
05
NCVXBQP1
10000
21
21
0.03
-
0.4922
E
+
08
0.5250
E
+
05
NCVXBQP2
10000
21
21
0.03
-
0.2812
E
+
08
0.3750
E
+
05
NCVXBQP3
10000
21
21
0.02
0.7034
E
+
07
0.3750
E
+
05
NOBNDTOR
5476
658
658
2.86
-
0.3683
E
+
14
0.3672
E
+
00
OBSTCLAE
10000
2780
2780
23.21
-
0.4897
E
+
12
0.1952
E
+
00
OBSTCLAL
10000
2274
2274
19.01
-
0.4304
E
+
12
0.1289
E
+
00
OBSTCLBL
10000
4281
4281
35.78
-
0.4744
E
+
12
0.1875
E
+
00
OBSTCLBM
10000
3998
3998
33.56
-
0.5849
E
+
12
0.6435
E
-
01
OBSTCLBU
10000
1551
1551
13.03
-
0.5328
E
+
11
0.4544
E
-
01
PENTDI
1000
14
14
0.01
-
0.4267
E
+
02
0.2967
E
-
05
TORSION1
10000
811
811
6.80
-
0.2733
E
+
14
0.2905
E
+
00
TORSION2
10000
1644
1644
13.76
-
0.5354
E
+
13
0.1543
E
+
00
TORSION3
10000
712
712
5.98
-
0.7571
E
+
11
0.4588
E
-
01
TORSION4
10000
2079
2079
17.47
-
0.3908
E
+
14
0.7154
E
+
00
TORSION5
10000
744
744
6.20
-
0.8637
E
+
13
0.1263
E
+
00
TORSION6
10000
730
730
6.07
-
0.1361
E
+
15
0.5556
E
+
00
TORSIONA
10000
831
831
7.65
-
0.3461
E
+
13
0.8291
E
-
01
TORSIONB
10000
1049
1049
9.61
-
0.4526
E
+
10
0.3910
E
-
01
TORSIONC
10000
1165
1165
10.65
-
0.6297
E
+
14
0.5303
E
+
00
TORSIOND
10000
980
980
9.01
-
0.3573
E
+
14
0.4092
E
+
00
TORSIONE
10000
704
704
6.48
-
0.2645
E
+
15
0.2775
E
+
01
TORSIONF
10000
646
646
5.91
-
0.1976
E
+
15
0.1459
E
+
01
ODNAMUR
11130
10692
10692
57.16
0.8793
E
+
04
0.1247
E
+
01
GRIDGENA
6218
1145
1145
6.92
0.2352
E
+
05
0.9992
E
-
03
NOBNDTOR
5746
—
—
—
—
—The symbol “—” indicates that the CUTEr system becomes nonrespondent when calculating the corresponding problem. For L-BFGS-B, we record the number of function and gradient evaluations Nfg because L-BFGS-B always evaluates the function and gradient at the same time.Some general observations on the results in Tables2 and 3 are the following. From the last column of Table 3, we see that, for 41 problems, L-BFGS-B has terminated abnormally without being able to satisfy the termination condition (32). SDPRP fails to reach a stationary point based on the stopping criteria (32) only on the 12 problems: BDEXP,EXPQUAD,SOND1LS,QR3DLS,BIGGS1,BQPGAUSS,CHENHARK,JNLBRNGB,NCVXBQP1,ODNAMUR,NCVXBQP2,GRIDGENA. It is an interesting fact that L-BFGS-B obtains at least as good function value as SDPRP on some problems but the projected gradient did not meet the stooping condition. This reason is not clear as pointed out by Zhu et al. in [21]. In analyzing the performance of both methods on the 13 problems at which at least one method works successfully, we observe that, on these problems, L-BFGS-B requires less iterations, less function evaluations than SDPRP. The sterical data of both algorithms are summarized in Table 4, where “Scueed” is the number of successes and “Failure” is the number of failures of tested problems for both algorithms.Table 4
Statical data of SDPRP and L-BFGS-B.
Success
Failure
SDPRP
42
12
L-BFGS-B
13
41Taking everything together, the preliminary numerical comparisons indicate that our proposed method is efficient and competitive with the well-known method L-BFGS-B. Moreover, we conclude that the method provides a valid approach for solving large-scale box-constrained problems.


## 5. Conclusions

In this paper, we have developed a subspace nonlinear conjugate gradient method for solving box-constrained optimization problems. For most of the tested optimization problems (42 out of 54), the algorithm works successfully to terminate at the solution. However, in the most cases, the number of function evaluations seems large. A common feather of this method is that all the generated points are feasible, without requiring gradient projection as usual methods in this scheme. Therefore, this may be the reason why the function evaluation numbers are higher than those of projected gradient or trust-region methods. However, we also believe that SDPRP is a valid approach for box-constrained problems. In our view, there are at least three issues that could lead to improvements. The first point that should be considered is probably the choice of the parameters in the active set identification technique and the value of the used parameters is not the only choice. Another important point that should be further investigated is the adoption of gradient projection technique. The third assumption is strong; how can we modify this algorithm so as to avoid the strict complementary assumption? Additionally, it is worthwhile to investigate that some existing conjugate gradient methods, such as [12, 22], whether it is possible to embed an active set framework to solve box-constrained problems. To this end, although the proposed method does not obtain significant development as we have expected, we think that the enhancement of this proposed method is still noticeable. Hence, we believe that the proposed algorithm is a valid approach for the problems and possibility and it may have its own potency.



---

*Source: 236158-2014-04-13.xml*