A Subspace, Interior, and Conjugate Gradient Method for Large-Scale Bound-Constrained Minimization Problems.
A subspace adaptation of the Coleman--Li trust region and interior method is proposed for solving large-scale bound-constrained minimization problems. This method can be implemented with either sparse Cholesky factorization or conjugate gradient computation. Under reasonable conditions the convergence properties of this subspace trust region method are as strong as those of its full-space version.Computational performance on various large test problems is reported; advantages of our approach are demonstrated. Our experience indicates that our proposed method represents an efficient way to solve large bound-constrained minimization problems.
Introduction
. Recently Coleman and Li [1, 2, 3] proposed two interior and reflective Newton
methods to solve the bound-constrained minimization problem, i.e.,
min
algorithms are interior methods since the iterates fx k g are in the strict interior of the feasible region, i.e.,
ug. These two methods differ in that a line search to update iterates is used
in [2, 3] while a trust region idea is used in [1]. However, in both cases convergence is accelerated with
the use of a novel reflection technique.
The line search method version appears to be computationally viable for large-scale quadratic problems
[3]. Our main objective here is to investigate solving large-scale bound-constrained nonlinear
minimization problems (1.1), using a large-scale adaptation of the Trust-region Interior Reflective (TIR)
approach proposed in [1].
The TIR method [1], outlined in FIG. 1, elegantly generalizes the trust region idea for unconstrained
minimization to bound-constrained nonlinear minimization. Here g k
. The crucial
role of the (diagonal) affine scaling matrices D k and C k will become clear in x2.
An attractive feature of the TIR method [1] is that the main computation per iteration is solving a
Research partially supported by the Applied Mathematical Sciences Research Program (KC-04-02) of the Office of Energy
Research of the U.S. Department of Energy under grant DE-FG02-90ER25013.A000, and by NSF, AFOSR, and ONR through
grant DMS-8920550, and by the Advanced Computing Research Institute, a unit of the Cornell Theory Center which receives
major funding from the National Science Foundation and IBM Corporation, with additional support from New York State and
members of its Corporate Research Institute.
y Computer Science Department, Cornell University, Ithaca, NY 14850.
z Computer Science Department and Center for Applied Mathematics, Cornell University, Ithaca NY 14850.
The TIR Method [1]
For
1. Compute define the quadratic model
2. Compute a step s k , with based on the subproblem:
3. Compute
4. If ae k ? - then set x
5. Update D k as specified below.
Updating Trust Region Size D k
1. If ae k - then set D k+1 2 (0;
2. If ae k 2 (-;
3. If ae k - j then
if
set
otherwise,
set
FIG. 1. The TIR Method for Minimization Subject to Bounds
standard unconstrained trust region subproblem:
min
The method of Mor- e and Sorensen [4] can be directly applied to (1.2) if Cholesky factorizations of
matrices with the structure of H k can be computed efficiently. However, this method is unsuitable for
large-scale problems if the Hessian H k is not explicitly available or (sparse) Cholesky factorizations are too
expensive. Recently, Sorensen [5] proposed a new method for solving the subproblem (1.2) using matrix
vector multiplications. Nonetheless, the effectiveness of this approach for large-scale minimization,
particularly in the context of our trust region algorithm, is yet to be investigated.
We take the view that solving the full space trust region subproblem (1.2) is too costly for a large-scale
problem. This view is shared by Steihaug [6] who proposes an approximate (conjugate gradient)
approach. Steihaug's approach to (1.2) seems viable although our computational experience (see Table
indicates that important negative curvature information can be missed, causing a significant increase
in the number of minimization iterations.
In this paper, we propose an alternative: an approximate subspace trust region approach (STIR). We
verify that, under reasonable conditions, the convergence properties of this STIR method are as strong
as those of its full-space version. We explore the use of sparse linear algebra techniques, i.e., sparse
factorization and preconditioned conjugate gradients, in the context of this approach.
In addition, we demonstrate the benefits of our affine scaling, reflection and subspace techniques
with computational results. First, for (1.1), our affine scaling technique outperforms the classical Dikin
scaling [7], at least in the context of our algorithm. Second, we examine our method with and without
reflection. We show the reflection technique can substantially reduce the number of minimization itera-
tions. Third, our computational experiments support the notion that the subspace trust region method is
a promising way to solve large-scale bound-constrained nonlinear minimization problems. Compared to
the Steihaug [6] approach, the subspace approach is more likely to capture negative curvature information
and consequently leads to better computational performance. Finally, our subspace method is competitive
with, and often superior to, the active set method in LANCELOT [8].
The paper is organized as follows. In x2, we briefly summarize the existing TIR method. Then we
provide a computational comparison of the subspace trust region method and the Steihaug algorithm in
the context of unconstrained minimization in x3. We introduce a subspace method STIR, and discuss its
convergence properties, in x4. Issues concerning the computation of negative curvature directions and
inexact Newton steps are discussed in x5; computational results are provided indicating that performance
is typically not impaired by using an inexact Newton step. Concluding remarks appear in x7. The
convergence analysis of the STIR method is included in the appendix.
2. The TIR Method. In this section we briefly review the full-space TIR method [1], sketched in
FIG. 1. This method closely resembles a typical trust region method for unconstrained minimization,
min x2! n f(x). The key difference is the presence of the affine scaling (diagonal) matrices D k and C k .
Next we briefly motivate these matrices and the TIR algorithm.
The trust region subproblem (1.2) and the affine scaling matrices D k and C k arise naturally from
examining the first-order Kuhn-Tucker conditions for (1): if a feasible point l ! x is a local minimizer,
then x i is not at any of its bounds. This characterization is
expressed in the nonlinear system of equations
where
and the vector v(x) 2 ! n is defined below: for each 1 - i - n,
(i). If
(ii). If
(iii). If
(iv). If g i - 0 and l
The nonlinear system (2.1) is not differentiable everywhere; nondifferentiability occurs when v
Hence we avoid such points by maintaining strict feasibility, i.e., restricting x k 2 int(F). A Newton step
for (2.1) is then defined and satisfies
where
Here J v n\Thetan corresponds to the Jacobian of jv(x)j. Each diagonal component of the diagonal
matrix J v equals to zero or \Sigma1. If all the components of l and u are finite, J
we define J v
Equation (2.3) suggests the use of the affine scaling transformation: -
x. This transformation
reduces the constrained problem (1.1) into an unconstrained problem: a local minimizer of (1.1) corresponds
to an unconstrained minimizer in the new coordinates -
x (for more details, see [1]). Therefore a
reasonable way to improve x k is to solve the trust region subproblem
min
where
s:
s. Subproblem (2.5) is equivalent to the following problem in the original variable space:
min
where
In addition to the close resemblance to an unconstrained trust region method, the TIR algorithm has
strong convergence properties with explicit conditions on steps for optimality. We now describe these
conditions.
The TIR algorithm requires strict feasibility, i.e., x use ff
to denote the
step obtained from d k with a possible step-back for strict feasibility. Let -
k denote the minimizer along
d k within the feasible trust region, i.e., -
ff
The above definition implies that '
Explicit conditions which yield first and second-order optimality are analogous to those of trust region
methods for unconstrained minimization [1]:
Assume that p k is a solution to min s2! nf/ k and fi q and fi q
are two
positive constants. Then s
Condition (AS.3) is necessary for first-order convergence; (AS.4), together with (AS.3), is necessary
for second-order convergence. Both conditions (AS.3) and (AS.4) are extensions of convergence conditions
for unconstrained trust region methods. In particular, when assumptions
are exactly what is required of trust region methods for unconstrained minimization problems.
Satisfaction of both conditions (AS.3) and (AS.4) is not difficult. For example, one can choose s k
so that / k (s k ) is the minimum of the values /
However, this does not lead to
an efficient computation process. In [3] and [2], we have utilized a reflection technique to permit further
possible reduction of the objective function along a reflection path on the boundary. We have found in [3]
and [2] that this reflection process significantly enhances performance for minimizing a general quadratic
function subject to simple bounds.
x
s
pr
FIG. 2. Reflection Technique
For all the computational results in this paper, s k is determined from the best of three points
corresponding to /
k [p R
k denotes the piecewise direction path with p k
reflected on the first boundary it encounters, see FIG. 2.
We can appreciate the convergence results for this approach by observing the role of the affine scaling
matrix D k . For the components x i which are approaching the "correct" bounds, the sequence of directions
f\GammaD \Gamma2
becomes increasingly tangential to these bounds. Hence, the bounds will not prevent a large
step size along f\GammaD \Gamma2
k g from being taken. For the components x i which are approaching the "incorrect"
bounds, f\GammaD \Gamma2
points away from these bounds in relatively large angles (the corresponding diagonal
components of D k are relatively large and g k points away from these bounds). Hence, a reduction of at
least
implies the scaled gradient fD \Gamma2
converges to zero (i.e., first-order optimality).
The scaling matrix used in our approach is related to, but different from, the scaling typically used
in affine scaling methods for linear programming. The affine scaling matrix D affine
commonly used in affine scaling methods for linear programming, is formed from
the distance of variables to their closest bounds. Our scaling matrix D 2
k equals to D affine
only when
j. (Note that even in this case we employ the square root of the quantities
used to define D affine
.)
Before we investigate a subspace adaptation of TIR, we demonstrate the effectiveness of our reflection
idea and affine scaling technique. We consider random problem instances of molecule minimization
[9, 10], which minimize a quartic subject to bounds on the variables. Table 1 and 2 list the average
number of iterations (over ten random test problem instances for each entry) required for the different
techniques under comparison. The notation ? in front of a number indicates that the average number is at
least this number because the iteration number exceeds 1000, the maximum allowed, for some instance.
The details of the algorithm implementation are given in x6.

Table

1 demonstrates the significant difference made by a single reflection. The only difference
With Reflection 34.1 41.7 66.8 83.4 93.6
Reflection 71.4 >210.1 >425.4 >302.2 > 408.5

The STIR algorithm with and without reflection: number of iterations
100 200 400 800 1000
unconstrained: D k 38.6 47.3 61.4 72.7 93.6
D affine
D affine
k >517.4 >617.6 >517.3 >1000 >1000

Comparison of the STIR scaling Dk and Dikin scaling D affine
number of iterations
between the rows with and without reflection is the following. Without reflection, s k is determined by the
best of the two points based on
determined by the best of
the three points based on /
k [p R
(with reflection). The superiority of using the
reflection technique is clearly demonstrated with this problem.
In

Table

2, we compare the computational advantage of the selection D k over D affine
the only
difference is the scaling matrix. We differentiate between problems that have an unconstrained solution
(no bounds active at a solution) and those with a constrained solution. We observe that, for unconstrained
problems, there is no significant difference between the two scaling matrices. However, for the constrained
problems we tested, the choice D k is clearly superior. We observe that when D k is used, the number
of iterations for a constrained problem is roughly the same as that for the corresponding unconstrained
problem. For D affine
k , on the other hand, the number of iterations for a constrained problem is much larger
than for the corresponding unconstrained problem.
3. Approximation to the Trust Region Solution in Unconstrained Minimization. There are two
possible ways to approximate a full-space trust region solution in unconstrained minimization.
Byrd, Schnabel, and Schultz [11] suggest substituting the full trust region subproblem in the unconstrained
setting by
min
where S k is a low-dimensional subspace. (Our implementation employs a two-dimensional choice for
Another possible consideration for the approximation of (1.2) is the Steihaug idea [6], also proposed
in the large-scale unconstrained minimization setting. In a nutshell, Steihaug proposes applying the
method of preconditioned conjugate gradients (PCG) to the current Newton system until either negative
curvature is revealed, the current approximate solution reaches the boundary of the trust region, or the
Newton system residual is sufficiently reduced.
We believe that a subspace trust region approach better captures the negative curvature information
compared to the Steihaug approach [6]. To justify this we have conducted a limited computational study
in the unconstrained minimization setting.
We implement the subspace method with the subspace S k defined by the gradient direction g k and the
output of a Modified Preconditioned Conjugate Gradient (MPCG) method applied to the linear Newton
. The output is either an inexact Newton step s IN
defined by,
or a direction of negative curvature, detected by MPCG. Algorithm MPCG is given in greater detail in
FIG. 11,

Appendix

B. Our implementation of the Steihaug method can also be found in Appendix B.
Both the Steihaug and subspace implementations are wrapped in a standard trust region framework for
the unconstrained minimization problem. For both methods the preconditioning matrix used is
where G k is the diagonal matrix computed from G k
The same strategy is used to update D k (see x6 for more details). We let D
used for the subspace method and k \Delta k G for the Steihaug method ([6]).
We used twenty different unconstrained nonlinear test problems. All but four are test problems
described in [12], but with all the bound constraints removed. The problems EROSENBROCK and
EPOWELL are taken from [13]. The last two problems, molecule problems MOLE1 and MOLE3, are
described in [9, 10]. For all problems, the number of variables n is 260. The minimization algorithm
terminates when kgk We use the parameter in both FIG. 11 and FIG. 12.

Tables

3 and 4 compare the Steihaug and subspace methods described above in terms of the number of
minimization iterations and the total number of conjugate gradient (CG) iterations. Table 3 shows problems
for which negative curvature was not detected, and Table 4 shows problems for which negative curvature
was detected. Although not included here, the function values and gradient norms (upon termination) were
virtually the same for both methods for all problems. Since these values were essentially the same among
the two methods, we only discuss the difference in iterations counts. The difference in minimization and
CG iteration counts is plotted in FIG. 3 and FIG. 4.
Most notable in Table 3 and the graphs of FIG. 3 is how strikingly similar the results are for the
Steihaug and subspace methods; the minimization with each method stops within two iterations of the
other in all cases. Furthermore, both methods take an identical number of total CG iterations except for
the problem BROWN1 where the Steihaug method takes four more iterations. When negative curvature
is encountered, shown in Table 4 and in FIG. 4, the iteration counts for each method are again similar
for a few problems. For most problems, however, the Steihaug method takes more iterations, and for
some problems the difference is substantial. This is particularly true for the problems CHAINWOOD,
MOLE1 and MOLE3 (for CHAINWOOD, problem 3 in FIG. 4, the total difference in iteration counts is
Minimization CG
Problem Subspace Steihaug Subspace Steihaug
1. BROWN1 27 29 39 43
2. BROWN3 6 6 6 6
3. BROYDEN1A 11 11 81 81
4. BROYDEN1B 5 5 34 34
5. BROYDEN2B 7 7 71 71
6. CHAINSING 22 22 188 188
7. CRAGGLEVY 21 21 125 125
8. DEGENSING 22 22 188 188
9. EPOWELL
10. GENSING 22 22 83 83
11. TOINTBROY 7 7 58 58
12. VAR 43 43 5590 5590

Comparison when only positive curvature is encountered: number of iterations
Minimization CG
Problem Subspace Steihaug Subspace Steihaug
1. AUGMLAGN 36 29 267 228
2. BROYDEN2A 22 19 247 196
3. CHAINWOOD 156 988 3905 3878
4. EROSENBROCK 44 46 52 86
5. GEROSE 23 33 166 165
6. GENWOOD 58 63 304 275
7. MOLE1 46 119 460 376
8. MOLE3 125 186 6311 5356

Comparison when negative curvature is encountered: number of iterations
positive curvature problems
minimization
iterations
excess Steihaug iterations
excess subspace iterations
positive curvature problems
iterations
FIG. 3. Comparison of subspace and Steihaug trust region methods for unconstrained problems
negative curvature problems
minimization
iterations
negative curvature problems
iterations
FIG. 4. Comparison of subspace and Steihaug trust region methods for unconstrained problems
explicitly noted as it is beyond the scale of the graph). In general the subspace method does take more CG
iterations on problems with negative curvature, but it is these extra relatively inexpensive CG iterations
that reduce the total number of minimization iterations. (Again, for the problem MOLE3 the difference
in CG iterations is explicitly noted in FIG. 4 as it is beyond the scale of the graph.)
A closer examination of the behavior of the two algorithms indeed shows that when negative curvature
is not encountered, both methods take similar steps. (In this case, if the trust region is large enough,
both methods in FIG. 11 and FIG. 12 will stop under the same conditions after the same number of
CG iterations, as displayed in Table 3.) By the nature of the algorithms, if the Steihaug method detects
negative curvature, then so will the subspace approach. However if the subspace algorithm detects
negative curvature, the Steihaug method may terminate before it finds negative curvature; and then it does
not converge (to a local minimizer) as quickly as the subspace method. The important role that negative
curvature plays is supported by the fact that the subspace method often moves in a substantial negative
curvature direction when the Steihaug method overlooks negative curvature. Furthermore, it is when the
trust region radius D k is small that the Steihaug method is most likely to stop early and miss negative
curvature. Thus it appears that the effectiveness of the Steihaug idea decreases as nonlinearity increases.
4. The STIR Method. Supported by the discussion in x3, we propose a large-scale subspace adaptation
of the TIR method [1] for the bound constrained problem (1.1).
In moving from the unconstrained subspace approach to the box-constrained setting, it seems natural
to replace the full trust region subproblem (1.2) by the following subspace subproblem
min
where S k is a small-dimensional subspace in ! n , e.g., a two-dimensional subspace. A two-dimensional
subspace for the trust region subproblem (2.5) can be selected from the span of the two vectors
k g and a negative curvature direction -
k . This suggests that we form S k from the
directions fD \Gamma2
g. Will such subspace formulations succeed in achieving optimality?
We examine this issue in more detail.
It is clear that the including the scaled gradient direction D \Gamma2
k in S k , and satisfying (AS.3), will
guarantee convergence to a point satisfying the first-order optimality conditions. Let us assume for now
that fx k g converges to a first-order point x   . To guarantee that x   is also a second-order point the
following conditions must be met.
Firstly, it is clear that a "sufficient negative curvature" condition must be carried over from the
unconstrained setting [14]. To this end, we can require that sufficient negative curvature of the matrix -
be captured if it is indefinite, i.e., S k must contain a direction w
w k such that
Secondly, it is important that a solution to (4.1) lead to a sufficiently large step - the potential
difficulty is running into a (bound) constraint immediately. This difficulty can be avoided if the stepsize
sequence, along the trust region solution direction, is bounded away from zero. Subsequently, we define:
DEFINITION 4.1. A direction sequence fs k g has large-step-size if lim inf k!1 jD 2
If fast local convergence is desired then the subspace S k should also contain a sufficiently accurate
approximation to the the Newton direction D \Gamma1
M k is positive definite and -
k . An
inexact Newton step -
k for problem (1.1) is defined as an approximate solution to
with accuracy
s IN
such that kr k k=k -
Can we select two-dimensional subspaces satisfying all three properties and thus guarantee quadratic
(superlinear) convergence to a second-order point? The answer, in theory, is yes - the subspace adaptation
of TIR algorithm (STIR) in FIG.5 is an example of a subspace method capable of achieving the desired
properties.
To ensure convergence to a solution, the solution sequence of the subspace trust region subproblems
(4.1) need to have large-step-size. Lemma 1 below indicates that this can be achieved if we set S
are two sequences of uniformly independent vectors in the sense
that lim inffkz each with large-step-size.
LEMMA 1. Assume that fw k g and fz k g have large-step-size with kD k w
Moreover, lim inf k!1 fkz the solution sequence fp k g to the subproblem (4.1) with
has large-step-size.
Proof. The proof is very straightforward and is omitted here.
For the STIR method, a natural extension of the condition (AS.4) necessary for second-order optimality
is the following.
Assume that p k is a solution to min s2! nf/ k and fi q and
are two positive constants. Then s
Theorem 2 below, with the proof provided in the Appendix, formalizes the convergence properties
of STIR.
THEOREM 2. Let the level set be compact and f
be twice continuously differentiable on L. Let fx k g be the sequence generated by the STIR algorithm in
FIG.5. Then
1. If (AS.3) is satisfied, then the Kuhn-Tucker condition is satisfied at every limit point.
The STIR Method
For
1. Compute define the quadratic model
2. Compute a step s k , with based on the subspace
subproblem,
where the subspace S k is set up as below.
3. Compute
4. If ae k ? - then set x
5. Update D k as specified in FIG.1.
Determine Subspace
[Assume that w
has large-step-size. Let
small positive constant.]
M k is positive definite
k is not positive definite
IF (D \Gamma2
(D \Gamma2
END
END
FIG. 5. The STIR Method for Minimization Subject to Bound Constraints
2. Assume that both (AS.3) and (AS.5) are satisfied and -
w k in FIG. 5 contains sufficient negative
curvature information whenever -
M k is indefinite, i.e.,
with
(a) If every limit point of fx k g is nondegenerate, then there is a limit point x   at which
both the first and second-order necessary conditions are satisfied.
(b) If x   is an isolated nondegenerate limit point, then both the first and second-order
necessary conditions are satisfied at x   .
(c) If -
is nonsingular for some limit point x   of fx k g and -
M k is
positive definite, then -
is positive definite, fx k g converges to x   , all iterations are
eventually successful, and fD k g is bounded away from zero.
The degeneracy definition is the same as in [1].
DEFINITION 4.2. A point x 2 F is nondegenerate if, for each index i:
We have established that in principle it is possible to replace the full-dimensional trust region
subproblem with a two-dimensional variation. However, the equally strong convergence properties of
STIR hinges on obtaining (guaranteed) sufficient negative curvature direction with large-step-size. We
discuss this next.
5. Computing Negative Curvature Directions with Large-Step-Size. Is it possible, in principle,
to satisfy both the sufficient negative curvature requirement (4.2) and the large-step-size property? The
answer is yes: let u k be a unit eigenvector of -
k corresponding to the most negative eigenvalue, i.e.,
. It is easily verified that for any convergent subsequence lim k!1 - min ( -
the sequence fD \Gamma1
has large-step-size.
However, it is not computationally feasible to compute the (exact) eigenvector u k . Therefore,
approximations, and short cuts, are in order. Can we compute approximate eigenvectors with large-step-
A good approximation to an eigenvector corresponding to an extreme eigenvalue can usually be
obtained through a Lanczos process [15]. Using the Lanczos method for -
k with an initial vector -
approximate eigenvectors at the j-th step are computed in the Krylov space
In the context of our algorithm, the vectors D \Gamma1
are natural choices for the initial vector
when applying the Lanczos method.
Our key observation is the following. If a sequence fD \Gamma1
k g has large-step-size then each sequence
in D \Gamma1
retains this property.
Now assume that -
w k is the computed vector from the Lanczos method which contains the sufficient
negative curvature information with respect to -
k . It can be verified, based on the recurrence relation,
that fD \Gamma1
k g all have large-step-size if the Lanczos vectors f -
retain orthogonality.
w k is in the Krylov space K( -
it is clear that fw
has large-step-size. In other
words, in order to generate a negative curvature direction sequence with large-step-size, orthogonality
needs to be maintained in the Lanczos process. Fortunately, as discussed in [16], it is quite reasonable
to assume that until all of the distinct eigenvalues of the original matrix have been approximated well,
orthogonality of the Lanczos vectors are well maintained. Since we are only interested in a direction with
sufficient negative curvature, we expect that it can be computed before loss of orthogonality occurs.
A second (and cheaper) strategy is to employ a modified preconditioned conjugate gradient scheme,
e.g., MPCG in FIG.12. Unfortunately, this process is not guaranteed to generate sufficient negative
curvature; nonetheless, as indicated in [17], the MPCG output will satisfy the large-step-size property.
Finally we consider a modified Cholesky factorization, e.g., [18], to obtain a negative curvature
direction. Assume that f -
is indefinite and fd k g is obtained from the modified Cholesky method. We
demonstrate below that fd
d k g has large-step-size under a nondegeneracy assumption.
The negative curvature direction -
computed from the modified Cholesky method (see
[18], page 111) satisfies
where L k is a lower triangular matrix, P k is a permutation matrix and e j k
is the j k th elementary vector,
(j is a bounded and non-negative diagonal matrix.
Without loss of generality, we assume that P I .
We argue, by contradiction, that fd k g has the large-step-size property. Assume that fd k g does not
have this property. From L T
and that L k is a lower triangular matrix with unit diagonals, it is
clear that -
Moreover, from -
definition (2.4) of -
k , the first components of fD k
are bounded. This implies that fv j k
converges to zero.
From the modified Cholesky factorization, the matrix -
is indefinite but -
is positive definite. But this is impossible for sufficiently large k because, again using the definition (2.4)
of -
converges to a matrix of the form
positive (because of the nondegeneracy assumption). Therefore, we conclude that fd k g has
large-step-size.
6. Computational Experience. We demonstrate the computational performance of our STIRmethod
given in FIG.5. Below we report our experience with the modified Cholesky and the conjugate gradient
implementations. We examine the sensitivity of the STIR method to a starting point. Finally,
some limited comparisons with SBMIN of LANCELOT [8] are also made.
In the implementation of STIR, we compute s k using a reflective technique as shown in FIG.2. The
exact trust region updating procedure is given below in FIG.6.
Updating Trust Region Size D k
1. If ae k - 0 then set D
2. If ae k 2 (0; -] then set D
3. If ae k 2 (-; j) then set D
4. If ae k - j then
if
otherwise
FIG. 6. Updating Trust Region Size
Our experiments were carried out on a Sun Sparc workstation using the Matlab environment.
The stopping criteria used are as follows. We stop if
either
or
or no negative curvature has been detected for -
We define
We also impose an upper bound of 600 on the number
of iterations.
We first report the results of the STIR method using the modified Cholesky factorization. Table 5
lists the number of iterations required for some standard testing problems (for details of these problems
see [12]). (For all the results in this paper, the number of iterations is the same as the number of objective
function evaluations.) The problem sizes vary from 100 to 10; 000. The results in Table 5 indicate that, for
these testing problems at least, the number of iterations increases only slightly, if at all, with the problem
size. Moreover, in comparison to the unconstrained problems, the presence of the bound restrictions does
not seem to increase the number of iterations. This is depicted pictorially in FIG. 7. In this graph, the
problem size is plotted versus iteration count. For each problem, the corresponding points have been
connected to show how the iteration count relates to the problem size.
Our second set of results are for the STIR algorithm but using a conjugate gradient implementation.
We use the algorithm MPCG in FIG.12 to find the directions needed to form the subspace S k . The
stopping condition applied to the relative residual in MPCG is 0:005. The results are shown in

Table

6 and FIG. 8. Again, for these problems the iteration counts are low and steady. The exception
is for the problem VAR C with 10; 000 variables, where the iteration count jumps to 86. This is one of
Problem 100 200 500 1000 10000
GENROSE C 11 11
GENSING U 24 25 25 26 27
GENSING C
DEGENSING C 28 28 28 28 29
GENWOOD
BROYDEN2A C 14 19 17 19 19
CRAGGLEVY U
VAR C

STIR method with exact Newton steps: number of iterations
problem size
iterations
FIG. 7. STIR performance with exact Newton steps
100 1000 1000050problem size
iterations
FIG. 8. STIR method with inexact Newton steps
several degenerate problems included in this test set. With a tighter bound j on the relative residual in
MPCG, we could decrease the number of minimization iterations for this problem (note that the STIR
with exact Newton steps only takes 38 iterations). However, this change would also increase the amount
of computation (conjugate gradient iterations).
Next we include some results which indicate that our STIR method is fairly insensitive to the starting
point. The results in Table 7 were obtained using exact Newton steps on problems of dimension 1000.
The results in Table 8 were obtained using the conjugate gradient implementation, also on problems with
1000 variables. The starting points are as follows: original is the suggested starting point according
to [12]; upper starts all variables at upper bounds; lower starts all variables at the lower bounds; middle
starts at the midpoint between bounds; zero starts each variable at zero (the origin); upper-lower starts the
odd variables at the upper and the even variables at the lower bounds; lower-upper is the reverse of this.
For all of these, we perturb the starting point slightly if necessary to be strictly feasible. Note that for the
problem BROWN3 C, the iteration count is not shown starting at middle and at origin as the gradient is
undefined at both these starting points. These results are also shown graphically in FIG. 9 and FIG. 10.
From these graphs it is clear that both implementations of STIR are fairly robust when it comes to starting
Problem 100 200 500 1000 10000
GENROSE U 21 21 21
GENSING U 23 23 24 24 25
GENSING C
CHAINSING U 21 21 21
CRAGGLEVY U 26 26
CRAGGLEVY C 26 26 26 26 27

STIR method with inexact Newton steps, krk=kgk - 0:005: number of iterations
Starting Point
Problem original upper lower middle zero up-low low-up
GENROSE C 11 27 33 15 16 43 27
CRAGGLEVY C 26 26 34 37

STIR method with exact Newton steps for number of iterations
points. This is in contrast to active set methods where the starting point can have a more dramatic effect
on the iteration count.
Last we contrast the performance of the STIR method using the conjugate gradient option with the
SBMIN algorithm, an active set method, in the LANCELOT software package [8]. In particular, we
choose problems where negative curvature is present or where it appears that the "active set" at the
solution may be difficult to find. We expect our STIR method to outperform an active set method in these
situations; indeed, we have found this to be the case. For these problems, we use the default settings for
LANCELOT and adjusted our STIR stopping conditions to be comparable if not more stringent.
First consider a constrained convex quadratic problem. The results, given in Table 9, show that our
proposed STIR method is markedly superior (by an order of magnitude) to SBMIN on this problem (c.g. it
is the total number of conjugate gradient iterations). SBMIN takes many iterations on this problem when
the starting point is near some of the bounds - the method mis-identifies the correct active set at the
solution and takes many iterations to recover. Our proposed STIR method, a strictly interior method,
moves directly to the solution without faltering when started at the same point.

Table

summarizes the performances of STIR and SBMIN, on a set of constrained problems
exhibiting negative curvature. (Again the problems are from [12] except the last two have been constrained
differently to display negative curvature.) STIR is significantly better on these problems - this is probably
due to the fact that negative curvature is better exploited in our subspace trust region approach than in
Starting Point
Problem original upper lower middle zero up-low low-up
GENSING C
28
DEGENSING C 33 43 37 42 37 37 44
GENWOOD C
28

STIR method with inexact Newton steps for number of iterations
orig up low mid zero up-lo lo-up50iterations
FIG. 9. STIR method with exact Newton steps at varied starting points
inexact STIR SBMIN
iteration c.g. it iteration c.g. it

STIR with inexact Newton steps vs. LANCELOT SBMIN on a convex quadratic: number of iterations
orig up low mid zero up-lo lo-up50iterations
FIG. 10. STIR method with inexact Newton steps at varied starting points
inexact STIR SBMIN
Problem 100 1000 10000 100 1000 10000
AUGMLAGN U 34
GENWOOD U 62 67 63 439 952 554
GENWOOD NC

STIR with inexact Newton steps vs. LANCELOT SBMIN when negative curvature exists: number of iterations
the Steihaug trust region method, which SBMIN employs. This is consistent with results presented in x3,
e.g., see

Table

4.
7. Conclusion. Based on the trust-region interior reflective (TIR) method in [1], we have proposed
a subspace TIR method (STIR) suitable for large-scale minimization with bound constraints on the
variables. In particular, we consider a two-dimensional STIR in which a subspace is formed from the
scaled gradient and (inexact or exact) Newton steps or a negative curvature direction.
We have designed and reported on a variety of computational experiments. The results strongly
support the different components of our approach: the "subspace idea", the use of our novel affine
scaling matrix, the modified Cholesky factorization and conjugate gradient variations, and the "reflection
technique". Moreover, preliminary experimental comparisons with code SBMIN, from LANCELOT [8],
indicate that our proposed STIR method can significantly outperform an active-set approach for some
large-scale problems.



--R

An interior
On the convergence of reflective Newton methods for large-scale nonlinear minimization subject to bounds
A reflective Newton method for minimizing a quadratic function subject to bounds on the variables.

Minimization of a large scale quadratic function subject to an ellipsoidal constraint.
The conjugate gradient methods and trust regions in large scale optimization.
Iterative solution of problems of linear and quadratic programming.
LANCELOT: A Fortran Package for Large-Scale Nonlinear Optimization (Release
The molecule problem: Determining conformation from pairwise distances.

A family of trust-region-based algorithms for unconstrained minimization with strong global convergence properties
Testing a class of methods for solving minimization problems with simple bounds on the variables.

Approximate solution of the trust region problem by minimization over two-dimensional subspaces
Matrix Computations.
Lanczos Algorithms for Large symmetric eigenvalue computations
Inexact Reflective Newton Methods for Large-Scale Optimization Subject to Bound Constraints
Practical Optimization.
--TR

--CTR
Detong Zhu, A new affine scaling interior point algorithm for nonlinear optimization subject to linear equality and inequality constraints, Journal of Computational and Applied Mathematics, v.161 n.1, p.1-25, 1 December
Jiaju Zheng , Shuying Cao , Hongli Wang , Wenmei Huang, Hybrid genetic algorithms for parameter identification of a hysteresis model of magnetostrictive actuators, Neurocomputing, v.70 n.4-6, p.749-761, January, 2007
L. N. Vicente, Local Convergence of the Affine-Scaling Interior-Point Algorithm for Nonlinear Programming, Computational Optimization and Applications, v.17 n.1, p.23-35, Oct. 2000
Amit Jain , David Blaauw, Slack borrowing in flip-flop based sequential circuits, Proceedings of the 15th ACM Great Lakes symposium on VLSI, April 17-19, 2005, Chicago, Illinois, USA
Manfred Weiler , Ralf Botchen , Simon Stegmaier , Thomas Ertl , Jingshu Huang , Yun Jang , David S. Ebert , Kelly P. Gaither, Hardware-Assisted Feature Analysis and Visualization of Procedurally Encoded Multifield Volumetric Data, IEEE Computer Graphics and Applications, v.25 n.5, p.72-81, September 2005
R. Deng , P. Davies , A. K. Bajaj, A nonlinear fractional derivative model for large uni-axial deformation behavior of polyurethane foam, Signal Processing, v.86 n.10, p.2728-2743, October 2006
interior method;trust region method;box constraints;bound-constrained problem;inexact Newton step;conjugate gradients;negative curvature direction