A Reflective Newton Method for Minimizing a Quadratic Function Subject to Bounds on some of the Variables.
We propose a new algorithm, a reflective Newton method, for the  minimization of a quadratic function of many variables subject to upper and lower bounds on some of the variables. The method applies to a general (indefinite) quadratic function for which a local minimizer subject to bounds is required and is particularly suitable for the large-scale problem. Our new method exhibits strong convergence properties and global and second-order convergence and appears to have significant practical potential. Strictly feasible points are generated. We provide experimental results on moderately large and sparse problems based on both sparse Cholesky and preconditioned conjugate gradient linear solvers.
Introduction
. In this paper we propose a new algorithm for solving the box-
constrained quadratic programming problem
The matrix H is symmetric and, in general, indefinite; l
We denote the feasible region and the strict
ug. When H is indefinite we are interested in locating
a local minimizer.
Problem (1.1) arises as a subproblem when minimizing general nonlinear functions
subject to bounds and as a problem in its own right. The box-constrained quadratic
programming problem represents an important class of optimization problems and has
been the subject of considerable recent work (e.g., [1, 5, 12, 13, 16, 19, 21, 24, 26, 27, 32]).
A special subclass deserves mention: the box-constrained least-squares problem,
where A is a rectangular m-by-n matrix with m ? n. Our proposed algorithm can
of course be applied to this special case if we form b. The
determination of a version of our algorithm which does not involve the formation of the
matrix H is an open question.
We propose a new approach, a reflective Newton algorithm. The algorithm generates
a sequence of strictly feasible iterates, fx k g, which converges under standard
assumptions to a local solution of (1.1), x   , at a quadratic convergence rate. Coleman
and Li [10] establish theoretical properties of the reflective Newton approach applied to
the general nonlinear box-constrained problem - as we indicate in Section 5 these results
apply directly to the reflective Newton procedure proposed here for the quadratic
minimization problem (1.1). In this paper we discuss the nature of the reflective transformation
(Section 2); we discuss the reflective Newton approach as applied to problem
with emphasis on a specialized line search to exploit the special structure of this
problem (Section 3); numerical experiments involving an implementation of the reflective
Newton method applied to box-constrained quadratic minimization (1.1) are
discussed (Section 4).
The sequence fx k g generated by the algorithm is strictly feasible: therefore, the
algorithm can be regarded as an "interior-point" algorithm. However, this is a misleading
classification. The algorithm differs markedly from methods commonly referred
to as "interior-point" algorithms. For example, the proposed algorithm does not use a
barrier function to ensure feasibility. The algorithm generates descent directions for q
and then follows a piecewise linear path, reflecting off constraints as they are encoun-
tered. Most interior point methods, on the other hand, generate descent directions (for
some function) and then restrict the step, along this straight-line direction, to ensure
feasibility.
The algorithm most similar to our current proposal is probably the recent method
due to Coleman and Hulbert [6]. (There is also a strong connection to previous work
by Coleman and Li [7, 8, 9, 20] on various norm minimization problems.) Both are
driven by the nonlinear system of equations representing first-order optimality condi-
tions. Both methods require piecewise quadratic minimization. The methods differ
in that our new algorithm is more general: positive definiteness of H is not required
and it is not necessary to have finite upper and lower bounds on all the variables -
the Coleman/Hulbert method requires both these restrictive properties. Finally, the
Coleman/Hulbert method is an exterior-point method, requiring strict decrease in a
piecewise quadratic "dual" function, whereas the new method generates feasible iter-
ates, requiring strict decrease in the original quadratic function q.
There are four key observations that underpin our new approach.
First, it is possible to change the constrained problem (1.1) to an unconstrained
problem without using a penalty parameter. We can replace (1.1) with an unconstrained
problem,
min
q(y)
where - q(y) is a continuous piecewise quadratic function of y, and y 2 R n is unrestricted.
Details of this transformation are given in the next section including a result, Theorem
1, proving the equivalence of (1.3) with (1.1). One view of our algorithm is that it is
designed to find a local minimizer of -
q(y). Variables x and y are related by a piecewise
linear transformation, a reflective transformation,
a feasible sequence fx k g. Moreover, evaluation of -
q(y) corresponds to evaluation of
q(R(y)).
Alternatively, one can view our approach entirely in the original variables x. Then,
instead of describing the method as a descent algorithm for the transformed problem
q(y), our method can be described as a method that generates feasible iterates by
following a piecewise linear path induced by the reflective mapping R. We discuss this
below.
The second key observation is that the first-order optimality conditions for (1.1),
or equivalently (1.3), can be expressed as a single system of nonlinear equations,
and a Newton step for this system is a descent direction for -
q in a neighbourhood of
a local solution y   . Moreover, in a neighbourhood of a local solution to (1.3) a full
Newton step for (1.4), i.e., a unit step size in the Newton direction, yields decrease in
q(y). This is a very important point because it suggests that a Newton process for (1.4)
is compatible with (1.3), at least in the neighbourhood of a solution. It suggests that
ultimate second-order convergence can be achieved while decreasing -
q(y).
The third observation leads to globalization of the Newton process. It turns out that
the Newton equation for (1.4), the nonlinear system representing first-order optimality
conditions, can be written in the form:
where M is a symmetric matrix. 3 Moreover, it turns out that M is positive definite in
a neighbourhood of a minimizer of -
and M can be interpreted, loosely, as a second
derivative matrix for -
q(y) . This suggests the use of an ellipsoidal constraint to ensure
a descent direction when far from the solution. Specifically, solve
\Deltag
where D is a positive diagonal scaling matrix and \Delta is positive. As we discuss in [10],
a good choice for matrix D(x) is
i.e., D is a diagonal matrix with the i th diagonal component equal to jv i (x)j 1
2 . Vector
defined in Figure 1 where
Diagonal matrix D plays an important role in this paper - henceforth we reserve the
notation D, without superscript, to refer to definition (1.7). Of course D k refers to (1.7)
with all quantities defined at the current point x k .
(iv) If g i - 0 and l
Fig. 1. Definition of v(x)
Solving (1.6) involves solving a symmetric positive definite system,
for a suitable -, and then s / D-s y . Thus it is easy to see that (1.6) leads to a descent
direction for - q at the current point. It may be felt that solving (1.6) is an expensive way
to determine a descent direction in the large-scale setting. With this is mind a restricted
version of (1.6) is used in our algorithm. In particular, similar to [2] we usually restrict
s to be in a low-dimensional subspace S. So (1.6) is replaced with
where S is a low-dimensional subspace of ! n . Provided the ellipsoidal constraint
inactive near the solution, and the low-dimensional subspace S
3 The function -
q is a piecewise quadratic function of y: therefore, r y -
q does not always exist. How-
ever, the proposed algorithm only generates points where r y -
q is defined.
ultimately includes the Newton direction, the solution to (1.9) will eventually be the
Newton step (1.5).
The fourth major ingredient of our approach is the line search. Once a descent
direction s k is determined, a one-dimensional line search is performed to approximately
locate a minimizer of -
q k has structure: -
q k is a one-dimensional
piecewise quadratic function and so an efficient specialized line search procedure can be
used. (Alternative view: A one-dimensional piecewise linear search is performed along
a "reflective path", p k (ff), defined by the reflective transformation R and beginning at
We conclude the introduction with a short review of optimality conditions for problem
(1.1).
The first-order optimality conditions can be written: If a feasible point x   is
a local minimizer of (1.1) then
Let F ree   denote the set of indices corresponding to "free" variables at point x
F ree
Second-order necessary conditions can be written 4 : If a feasible point x   is a local
minimizer of (1.1) then D 2
ree   is the submatrix of H
corresponding to the index set F ree
These conditions are necessary but not sufficient. To state practical sufficiency
conditions we first need a definition of degeneracy.
Definition 1. A point x nondegenerate if, for each index i:
With this definition we can state second-order sufficiency conditions: If a
nondegenerate feasible point x   satisfies D 2
ree   ? 0, then x   is a local
minimizer of (1.1).
2. The Reflective Transformation. One interpretation of our approach to solving
the box-constrained quadratic programming problem (1.1) involves a transformation
to an unconstrained piecewise quadratic minimization problem (1.3). The purpose of
this section is to introduce this transformation. Since some of the ideas involved are
more generally applicable, we begin our discussion at a more abstract level and gradually
work our way back to the box-constrained quadratic programming situation.
4 Notation: If a matrix A is a symmetric matrix then we write A ? 0 to mean A is positive definite;
means A is positive semi-definite.
Consider the problem
where f is a continuous and C is a closed connected region of ! n .
We consider when the constrained problem (2.1) can be replaced with an unconstrained
minimization problem of the form,
min
where R is a continuous onto mapping,
What further restrictions on the mapping R make this an acceptable transfor-
mation? To see that continuity is not enough consider the following one-dimensional
example. Let Obviously there is only one local solution
(the global solution), x   = 1. However, let R(y) be any continuous function, mapping
onto [0; 1], with a strict local maximizer at -
y with -
1). It is easy to
see that -
y is a local minimizer of f(R(y)) but - x is clearly not a local solution to the
original problem.
The following property plays the key role in answering this question.
R is an open mapping if for each ffl ? 0 and pair
See Munkres [25], for example, for a discussion of open mappings. We can now
answer our question.
Theorem 1. C be a continuous onto mapping. Further, assume R
is an open mapping. Then,
(i) If y   is a local minimizer of (2.2) then -
local minimizer of (2.1).
(ii) If x   is a local minimizer of (2.1) then for each -
- y is a local minimizer of (2.2). Moreover, there exists at least one -
y such that
(iii) Problem (2.1) is unbounded below if and only if problem (2.2) is unbounded
below.
Proof. (i) Assume y   is a local minimizer of (2.2). Let -
local minimizer, there exists ffl ? 0 such that
But by (2.3) there exists
Hence for each x   2 N
Therefore, - x is a local minimizer of (2.1).
(ii) Assume x   is a local minimizer of (2.1). But R is an onto mapping and therefore
there exists - y is a local minimizer, there exists ffl ? 0
such that
By continuity there exists
Therefore, for all y
Hence, - y is a local minimizer of (2.2).
(iii) Suppose fy k g is a sequence such that
lim
Alternatively, assume fx k g 2 C and
lim
But R is an onto mapping; therefore, for each x k there exists y k such that R(y k
and
lim
and therefore (iii) is established.
To illustrate, consider the problem
0g. A definition of R that clearly satisfies the open mapping property
is multiplication. Note that R is
differentiable and the Jacobian of R, J R (y), is nonsingular if and only if R(y) 2 int(C),
where int(C) is the interior of C. Specifically, r y
x
where D g is the diagonal matrix diag(r x f) and D y is the diagonal matrix diag(y). (Note
that r 2
y R is a tensor term and r x fr 2
y R is a matrix - diagonal in this case.) Therefore,
this definition of R leads to an unconstrained twice-differentiable minimization problem
and standard techniques can be used to solve (2.2). Unfortunately, our numerical
experience with this approach has been mixed: In particular, as problems become large
and ill-conditioning (and near-degeneracy) increases, the number of iterations required
by standard minimization algorithms, to achieve good accuracy, becomes quite large.
We feel this is due in part to the fact that this transformation causes an increase in
the complexity of the function to be minimized: e.g., a quadratic function becomes a
quartic. Our objection to this approach is largely numerical - ill-conditioning in the
original problem is accentuated when the problem is transformed to a more complex
form. Note also that the transformed problem may have many more local minimizers
- by Theorem 1 this, in itself, is not a problem. However, along with this increase in
the number of local minimizers comes an increase in negative curvature and this may
cause some optimization algorithms some difficulty. In any event, our experience with
this simple differentiable transformation has not been satisfactory: the subject of this
paper is an alternative definition of R.
For problem (2.5) consider is a vector, jvj denotes the vector
whose components are the absolute values of the vector v. It is clear that the open
mapping property holds. Note that R is not everywhere differentiable. In particular,
R is differentiable at point y if and only if R(y) 2 int(C), i.e., y i 6= 0. In this case
obviously nonsingular. Using this transformation,
f(R(y)) has a piecewise differentiable nature as a function of y. For example, if f(x) is
a quadratic function then f(R(y)) is a piecewise quadratic function of y.
We now extend the absolute value approach, to handle the more general
situation,
where for each index i either u i is finite or u Similarly, for each index i either l i
is finite or l
For simplicity we assume that the finite values of u are all equal to unity and the finite
values of l are all equal to zero (a simple translation and scaling can achieve this form 5 ).
The transformation we propose, x = R(y), is a diagonal transformation, i.e., for
5 A definition of the reflective transformation applied directly to the general problem is given in [10]
y
Fig. 2. The 1-Dimensional Reflective Transformation (Finite Upper and Lower Bound)
each index i, x i depends only on y i . This transformation, induces a piecewise
linear "reflective" path in x.
For example, if u Figure 2.
If l is the absolute value function; if l
then R i is illustrated in Figure 3.
Finally, if u . The four cases are described
more formally in Figure 4.
It is easy to verify that R satisfies the requirements of Theorem 1; therefore, use
of R does not introduce extraneous local minimizers nor does it restrict the set of local
minimizers.
Using the reflective transformation, problem (1.1) can be replaced with the unconstrained
piecewise quadratic minimization (1.3). In principle, problem (1.3) can be
solved by a descent direction algorithm, e.g., Algorithm 1 in Figure 5.
An advantage of using this piecewise linear transformation R is the linear aspect
of the transformation: when y is a differentiable point the local complexity of
is the same as the local complexity of f(x). When
q is a piecewise quadratic function. The apparent disadvantage is the
piecewise nature of -
f(y). This lack of differentiability means that conventional nonlinear
minimization methods cannot be used.
In particular, in order to guarantee convergence, restrictions on the nature of the
descent direction s y must be imposed. To see this suppose that y k is very close to a
y is a descent direction for - q at y k . If s y is nearly
perpendicular to this hyperplane then the usual descent condition, r y -
only result in a very small step since r y -
q is not continuous at x In [10]
we describe two properties, "constraint-compatibility" and "consistency", which help
guarantee sufficiently long steps and, consequently, global convergence. We discuss this
y
OE
x
Fig. 3. The 1-Dimensional Reflective Transformation with Infinite Upper Bound
briefly in Section 3.
The straight-line direction s y
k corresponds to a piecewise linear path in x. This
piecewise linear path can be described as follows. For simplicity, and without loss of
generality, assume y
k . Define the vector 6
Component i of vector BR k records the positive distance form x k to the breakpoint
corresponding to variable x k i in the direction s x
k . The piecewise linear (reflective) path
is defined by Algorithm 2 in Figure 6. Since only a single outer iteration is considered,
we do not include the subscript k with the variables in our description of Algorithm 2
- dependence on k is assumed.
Given the current point x k and a descent direction s x
denote the piecewise
linear path defined by Algorithm 2: For fi
Note that it is now possible to describe Algorithm 1 entirely in x-space without
explicitly introducing either the function -
q or the variables y. We do this in Algorithm

Figure

8.
The difference between Algorithm 1 and Algorithm 3 is notational. The view
presented by Algorithm 3 has the advantage that it is in the original space - visualization
6 For the purpose of computing BR we assume the following rules regarding arithmetic with infinities.
If a is a finite scalar then a
Case 1: (l
To evaluate x
then we can differentiate R to obtain the i th diagonal element of
the diagonal Jacobian matrix
else J R
Case 2: (l
To evaluate x
then we can differentiate R to obtain the i th diagonal element of the
Jacobian matrix
Case 3: (l
To evaluate x
then we can differentiate R to obtain the i th diagonal element of the
Jacobian matrix
If y
Case 4: (l
In this case there are no constraints on x i and so x
Fig. 4. The Reflective Transformation R
Algorithm 1
Choose
For
1. Determine a descent direction s y
for -
q(y) at y k
2. Perform an approximate line minimization of -
k ), with respect to
ff, to determine an acceptable stepsize ff k (such that ff k does not correspond
to a breakpoint)
3. y
Fig. 5. Descent dir'n algorithm for -
f (y)
Algorithm 2 [Let fi
[i u is a finite upper bound on the number of segments of the path to be determined]
For
1. Let fi i be the distance to the nearest breakpoint along
2.
3. Reflect to get new dir'n and update BR:
(a)
(b) For each j such that (b i
Fig. 6. Determine the linear reflective path p
of the reflective process is natural. The advantage of the first view, Algorithm 1, is that
the algorithm is a straight line descent direction algorithm, a familiar structure. It is
probably useful for the reader to keep both views in mind. In this paper we will primarily
work in the x-space and Algorithm 3. For simplicity we now drop the superscripts y
and x (e.g., s x becomes s).
It is well known that a descent direction algorithm demands sufficient decrease at
every step in order to achieve reasonable convergence properties. We use conditions
suggested by Goldfarb [18] for use in the unconstrained setting: Given
and a descent direction s k , ff k satisfies our approximate line search conditions if
and

Fig. 7. A Reflective Path
Algorithm 3
Choose
For
1. Determine an initial descent dir'n s x
k for q at x k 2 int(F ). Determine the
piecewise linear reflective path p k (ff) via Algorithm 2.
2. Perform an approximate piecewise line minimization of q(x k +p k (ff)), with
respect to ff, to determine an acceptable stepsize ff k (such that ff k does
not correspond to a breakpoint).
3. x
Fig. 8. A reflective path algorithm
is the
piecewise linear path defined by (2.8). Condition (2.9) can be interpreted as restricting
the step length from being too large relative to the decrease in f ; condition (2.10) can
be interpreted as restricting the step length from being relatively too small.
A basic reflective path algorithm for problem (1.1) can now be stated, Algorithm
4. To allow for flexibility, especially with regard to the Newton step, we do not always
require that both (2.9) and (2.10) be satisfied. Instead, we demand that either both
these conditions are satisfied or (2.9) is satisfied and ff k is bounded away from zero,
e.g., latter conditions are used to allow for the liberal use of Newton
steps and do not weaken the global convergence results.
Note that Algorithm 4 generates strictly feasible points; i.e., since x 1 2 int(F ), it
follows that x k 2 int(F ).
Algorithm 4 [ ae is a positive scalar.]
Choose
For
1. Determine an initial descent dir'n s k for q at x k . Note that the piecewise
linear path p k is defined by x
2. Perform an approximate piecewise line minimization of q(x k +p k (ff)), with
respect to ff, to determine ff k such that:
(a) ff k does not correspond to a breakpoint
(b) condition (2.9) is satisfied
(c) Either
i. ff k satisfies condition (2.10), or
ii.
3. x
Fig. 9. A reflective path algorithm satisfying line search conditions
3. Algorithm Specifics. A framework for our reflective Newton approach was
presented in the previous section, Algorithm 4. In this section we specify more precisely
how the search directions will be generated as well as the mechanics of the line search,
specialized to the quadratic problem (1.1).
The convergence analysis given in [10] uses two important properties of the sequence
of search directions, "constraint-compatibility" and "consistency". "Constraint-
compatability" is needed to guarantee that a sufficiently long step is taken before the
first constraint is encountered. The usual descent condition, that g T
sufficiently
negative, is not enough in the context of a reflective algorithm because this condition
takes no account of the proximity of the constraints. "Consistency" is a more standard
notion capturing the idea that first-order descent, represented by the term g T
sisitent with first-order convergence. Following (1.7), define D 2
Definition 2. A sequence of vectors fw k g is constraint-compatible if the sequence
fD \Gamma2
is bounded.
Definition 3. A sequence of vectors fw k g satisfies the consistency condition
Central to our approach, both in terms of achieving quadratic convergence and the
satisfaction of constraint-compatibility and consistency, is the frequent use of a reduced
trust region problem to determine s
where S k is a subspace of R n , D k is a positive diagonal scaling matrix, and
The matrix M k is defined:
where J v is the Jacobian 7 of v, where v is defined in Figure 9. Matrix D g
v is a diagonal
matrix with component i defined D
is an "extended
gradient", extended to deal with possible degeneracy. In particular,
where - g is a small positive constant. Clearly if x is a nondegenerate point and - g is
sufficiently small then (which implies that x is
degenerate) then
The diagonal matrix D(x), used in (3.1), is defined by (1.7), i.e. 8 ,
This choice yields a well-defined trust region problem (3.1). To see this note that using
(3.4), (3.1) becomes
where
s;
and
is a diagonal matrix, D
M k is positive definite in a
neighbourhood of a nondegenerate point satisfying the second-order sufficiency condi-
tions. Moreover, unlike fM k g, f -
k g is bounded. Matrix -
M k is a featured performer in
our reflective Newton algorithm. A Newton step is defined when -
M k is positive definite:
A final remark on the choice of scaling matrix (3.4). If we assume that D has the
1is the only reasonable choice. To see this suppose
consider that
7 Matrix J v is a diagonal matrix with each diagonal component equal to zero or unity. For example,
if all the components of u and v are finite then J has a finite lower bound and an
infinite upper bound (or vice-versa) then strictly speaking v i is not differentiable at a point
define J v
at such a point. Note that v i is discontinuous at such a point but v i \Delta g i is continuous.
8 Notation: If z is a vector then jzj 1
2 denotes a vector with the i th component equal to jz
and the calculation of DMD involves division by jv(x)j 1\Gamma2p
which includes components which go to zero as x ! x   . On the other hand, if approaches singularity as x ! x   (consider v
We will specify subspace S k below; it is important to realize that the cardinality of
our implementation. Therefore, the cost of solving (3.1)
is negligible. Given S k , the subspace trust region problem (3.5) can be approached
in the following way. Let S k be defined by the t k independent columns of an n-by-t k
be an
orthonormalization of the columns of D \Gamma1
for some vector s Y k . Therefore problem (3.5) becomes
ks Y k
and set s
. The solution to (3.7) is of negligible cost once the matrices are
small (see Appendix).
Algorithm 5 in Figure 10 presents a second-order reflective path algorithm.
Algorithm 5
Choose
For
1.
2. Determine initial descent dir'n s k for q at
M k is positive definite
and k-s N
k . If -
k is not positive definite choose
choose subspace S k , and solve (3.1) to get s k .
3. Determine ff
k and x k
otherwise, perform an approximate piecewise line minimization of q(x k
with respect to ff, to determine ff k such that
(a) ff k is not a breakpoint
(b) ff k satisfies (2.9) and (2.10).
4. x
Fig. 10. A second-order reflective path algorithm
Note: If ff accepted by the line search but corresponds to a breakpoint, then
modify
where ~
ff k is not a breakpoint, ~
9 If A is a matrix then ! A ? denotes the space spanned by the columns of A.
It remains to be more precise about the determination of s k and S k and to fully
specify the line search. We begin with s k and S k . Algorithm 6 in Figure 11 describes
our procedure.
Algorithm 6 [Let positive constants.]
Case 0: -
M k is positive definite and k-s N
Case 1: -
M k is positive definite and k-s N
if kr(-s N
solve (3.1) to get s k .
else
set s
Case 2: -
k is not positive definite. Compute
w k is a unit vector
such that fw k g is constraint-compatible and
z
solve (3.1) to get s k .
else
solve (3.1) to get s k .
Fig. 11. Determination of the descent direction s k
Remark on Case 2: We determine an appropriate negative curvature direction in
the following way. If a (sparse) Cholesky factorization of -
k does not complete then
k is not positive definite and a unit direction of non-positive curvature, -
readily
available and easily computable (e.g., [17]). Algorithm 6 can make use of -
sufficient negative curvature is displayed by -
and fw
is constraint-compatible. A constraint-compatibility test is implemented
by introducing a large constant, - cp , and requiring,
If either condition (3.9) or condition (3.10) is not satisfied then -
must be rejected. In
this case we can turn to a Lanczos process [10] to get a unit vector -
w k such that both
and (3.10) are satisfied. It is interesting to note that in our extensive numerical
experimentation, with results reported in Section 4, conditions (3.9) and (3.10) were
always satisfied by the (partial) Cholesky factorization method - the backup Lanczos
procedure was never required.
The Line Search. We have designed a specialized approximate line search procedure
to efficiently exploit the structure of this problem and to guarantee the line search
conditions in Algorithm 5. Before describing the approximate procedure, we develop
an exact line search procedure - this is possible because the problem is to find a local
minimizer of a quadratic function along a piecewise linear path. In the end we do not
use the exact line search per se but rather we use a truncated version of it, subroutine
"improve", within an overall approximate strategy. But we begin with the exact search.
The Exact Line Search. We are initially concerned with the exact determination
of ff
k where ff
k is a local minimizer of q(x k (ff)). Note: It is convenient to describe
the exact line search in terms of the y-variables, i.e., y
R(y view we have a straight-line minimization of a
piecewise quadratic function -
Henceforth in this section we omit the major iteration subscript.
The function -
is a continuous piecewise quadratic function. The ray
can be divided into intervals from left to right, I
q(ff) is smooth on each interval. Denote the restriction
of -
q(ff) to the j th interval by q j (ff) - note that q j (ff) is a quadratic function of a single
variable.
Our exact line search algorithm visits the intervals I 1 ; I 2 ; :::; in a left-to-right fashion
in an attempt to locate the first local minimizer of -
Assume we have not located
a local minimizer on intervals I 1 ; I 2 ; :::; I j \Gamma1 and assume that (q j There are
two possibilities: either q j (ff) has a minimizer strictly within the interval I j or it does
not. If it does, i.e., ff j
. However, if
does not admit a minimizer within I j then we must consider the possibility that fi j
is a minimizer of -
(ff). This is now the case if (q j+1
is not in int(I j )
and (q j+1 process is repeated on interval I j+1 .
Algorithm 7 in Figure 12 presents a compact description of the exact line search
algorithm we have sketched above.
Step (1.2) in Algorithm 7 follows from the observation that if (ff k
then s is a direction of infinite descent for -
q, beginning at y and therefore, by Theorem
descent for (1.1).
Algorithm 7 [Exact line search along direction s beginning at point y]
(0) Determine the array of breakpoints BR according to (2.7). fi 0 / 0.
(1) For
, the minimizer of q k , if it exists; otherwise, set ff k
1.
(problem (1.1) is unbounded).
, exit.
exit.
is the index such that according to
Algorithm 2.
Fig. 12. The Exact Line Search Algorithm
Our final concern, with regard to the exact line search, is an efficient implementation
of step (1.1). In theory this computation is straightforward. Assume q k
a k
. If a k
unbounded below; if a k
then q k (ff) is unbounded below (unless a k
which case q k is constant). However,
the challenge is to determine a k
1 and a k
efficiently, for (Note that a k
0 is not
The key to efficiency here is the observation that the reflective transformation R
is linear on each interval I is constant within each interval. For
each interval I k , define oe k to be the vector of diagonal elements of J R evaluated at any
point in the interval (y it follows that
for It follows that
Therefore, in the terminology used above, a k
a k
A straightforward implementation for determining a k
breakpoint (a k
0 is not needed). However, there is considerable structure that can be
exploited; in particular, oe k+1 , a vector with each component equal to \Sigma1, differs from
oe k in exactly one component. We can exploit this to reduce the work in the line search
to O(n) per breakpoint.
Suppose we have determined that I k does not contain ff k
and we have at hand the
following quantities:
a k
a k
and x k , the value of x at the k th breakpoint. If j is the index such that
then
The vector w k can be updated as follows:
where H j is column j of H. Coefficient a k+1
2 is simply computed:
a k+1
Coefficient a k+1
1 can be efficiently computed by considering the following equalities:
a k+1
Finally, x k+1 is computed:
In summary, the coefficients a k+1
2 and the intermediate quantities x
can be computed, given a k
using (3.16),(3.17),(3.18) and (3.19). This amounts
to approximately 4n work. Of course the initial quantities, w
a 0must be computed
from scratch requiring O(n 2 ) work. Therefore, if k br denotes the number of breakpoints
crossed, the total cost of the exact line search is:
initialization of w
(ii) O(n) for determination of BR,
(iii) O(k br n) for steps (1.0) - (1.5).
The Approximate Line Search. The exact line search described above is not
practical, for two reasons. First, an exact minimizer along a line may correspond exactly
to a breakpoint, i.e., a boundary point, and the algorithm requires strictly feasible
points. This is actually not a serious problem since a small perturbation would yield
strict feasibility and the reflective Newton method is not very sensitive to boundary
proximity.
A more serious objection to the exact line search is economy: despite the efficient
implementation described in the previous section, the relative cost can be high if the
number of breakpoints crossed, k br , is large. Certainly if there are a large number of
tight variables at the solution, say something close to n, then the total cost of the exact
line search algorithm ultimately becomes O(n 2 ) per line search. This is unsatisfactory
and unnecessary since an economical approximate line search can be just as effective.
In this section we describe an approximate line search, henceforth refer to as subroutine
improve, which uses the exact line search, described above, in a limited fashion-
beginning at an approximate minimizer, subject to a bound, k u , on the number of
breakpoints permitted to cross. In particular, improve is used in a cleanup role: after
determining an initial approximate minimizer by a bisection strategy, improve is called
upon to apply the exact line search strategy. Thereby the approximate minimizer is
further improved at cost O(k u n), where k u is typically chosen to be small, e.g., k
(In improve we also impose an approximate upper bound ff max on the size of ff. That
is, the size of the improvement is bounded by ff
Subroutine improve has the following calling sequence:
A more precise description of the approximate line search algorithm is given in

Figure

13, Algorithm 8. The basic idea is as follows. First, if the direction s k is a
Newton direction s N
k then a unit step is attempted. If (2.9) is satisfied then the full
Newton step is accepted subject to further improvement by subroutine improve and
possible (slight) adjustment to avoid a breakpoint. Second, if s k does not corrspond to
a Newton direction or if it does but a unit step does not satisfy (2.9), then a bisection
procedure is executed on the interval (0; ff u ) where ff u - 1 is an upper bound on the
step size. A point is located satisfying both (2.9) and (2.10) and then possibly further
improved with subroutine improve.
Algorithm 8 [Approximate line search along direction s beginning at point y]
k and a unit step along s k satisfies (2.9)
ff k corresponds to a breakpoint, set ff
is not a breakpoint
else f s k 6= s N
k or a unit step along s N
k does not satisfy (2.9)g
ffl use bisection to find -
and (2.10) such that -
ff k is not
a breakpoint
set ff
ff k corresponds to a breakpoint, set ff
is not a breakpoint
Fig. 13. Approximate Line Search Algorithm
4. Numerical Experiments. We have implemented our algorithm in a version
of Matlab which allows for sparse matrix data structures [15], now Matlab 4:0. In this
section we present some preliminary numerical results.
With the exception of the results reported on Table 12, all experiments were performed
on Sun Sparc workstations in the Matlab environment [22]. Experiments reported
in Table 12 were performed in a heterogeneous environment involving an Intel
IPSC/860 32-node multiprocessor as the "backend", and a Sun Sparc workstation as the
"frontend". Matrix factorizations and solves were performed on the "backend", in C,
while the main Matlab program executed on the "frontend". Communication between
"frontend" and "frontend" over ethernet was implemented through the use of Matlab
"MEX" files. We used this environment to facilitate the solving of very large problems.
(Details on this heterogeneous environment are given in [3].)
Starting and Stopping: In all the experiments reported in this paper the starting value
of x, i.e., x 1 , is as follows. For component j where both upper and lower bounds are
finite, choose the midpoint, j. If both upper and lower bound corresponding
to component j are infinite in size, choose choose
(Note: The
reflective Newton approach is not particularly sensitive to starting value. For example,
we repeated many of the experiments reported here using a random (strictly) feasible
starting point - very little difference in behaviour was detected.)
Choosing a robust stopping rule in optimization is not easy. Our primary stopping
rule is based on the relative difference in function value. This is reasonable partly
because strict feasibility is always maintained, and partly because often the real objective
in practical optimization is to achieve a point of relatively low function value.
Specifically, are primary stopping rule is:
We choose in Matlab on
a Sun Sparc workstation, . We do have secondary stopping criteria
as well - designed to determine when progress is deemed too slow. This secondary rule
tends to kick in when solving degenerate or ill-conditioned problems and a very flat
region around the solution has been entered.
Parameter settings: There are a number of parameters in the algorithm: most are
either in the very large or very small category. Here are the settings we used in our
experiments:
Used in the determination of scaling matrix D, see (3.3): -
Used in the line search, see (2.9): oe l = :1
Used in the line search, see (2.10): oe
A bound on the number of breakpoints crossed in subroutine improve:
ffl ae: A lower bound on the stepsize, see Algorithm 8:
If the line search produces a unit step which turns out to be a breakpoint,
this point is perturbed by an amount bounded by - ff kD k g k k, see (3.8): -
Used to test for constraint-compatibility, see (3.10): -
Used in the negative curvature test, see (3.9): ffl
Used in Algorithm
Used in Algorithm 8:
An upper bound on the bisection process used in Algorithm 8: ff
4.1. Positive Definite Problems. We have generated a number of quadratic
test problems with certain properties. In the first set of results we concentrate on
the case where H is symmetric positive definite. In the results reported below we use
sparse matrices H with sparsity patterns representing 3-dimensional grid using a 7-
point difference scheme. The Mor'e/Toraldo [24] QP-generator was adapted to generate
problems with a given sparsity pattern (see also [6]). We will not review the generator
characteristics here: our generator is a straightforward adaptation of the Mor'e/Toraldo
scheme to the sparse setting. We use several sparse Matlab functions (e.g., "sprandsym",
"sprand").
In

Tables

1-3, the dimension of the test problems is in each case. The
parameter "pctbnd" indicates the percentage of variables tight at the solution - approximately
evenly divided between upper and lower bounds. Parameter "deg" reflects
the degree to which the solution is (nearly) degenerate - the larger the value of "deg",
the greater the amount of (near) degeneracy. Technical details of "deg" are discussed in
[6]. Parameter "cond" reflects the conditioning of the matrix H: the condition number
of H is approximately 10 cond .
The upper and lower bound vectors, u and l, were generated as follows. Approximately
75% of the components of l were chosen to be finite and assigned the value of zero
- the index assignment was made in a random fashion. Similarly, approximately 75%
of the components of u were chosen to be finite and assigned the value of unity. Again,
the index assignment was made in a random fashion, independent of the assignment of
l.
Each row of Tables 1-3 reflects the results of 10 independent runs with the same
parameter settings. The third column, labelled "max", indicates the maximum number
of iterations required, over the set of 10 independent runs, to achieve the stopping
criteria; the fourth column, labelled "avg" records the average number of iterations
required to reach the stopping criteria over the 10 problems; the last column, labelled
"acc", records the number of digits of accuracy achieved in the function value (the true
solution is known).

Table
Positive
deg cond max avg acc
9 6 15 15 15
9 9
Observations on Tables 1-3: First, we observe the remarkable consistency of our reflective
Newton method on these problems. In terms of iterations required to achieve the
stopping criteria and accuracy attained in the function value, there is apparently very
little sensitivity to degeneracy, conditioning, or number of variables tight at the solution.
Of course we do not claim that accuracy in x is independent of condition/degeneracy
it surely is not. However, it is usually acceptable in optimization to attain a point
with nearly optimal function value and we have been quite successful in that (on this
test collection).
Second, the absolute number of iterations required to obtain a very accurate solution
(in terms of the function value q) is modest in every case, i.e., less than 20.
Positive
deg cond max avg acc
9 6 17 17 15
6 9 17 17 15
9 9 17 16.3 15

Table
Positive
deg cond max avg acc
9 3 17 16.3 15
9 6
Positive Definite Problems: Timing Breakdown
1000
This is very encouraging considering the dimension of the problems (n = 1000) and the
spectrum of problem characteristics being considered.
It is important to know where the algorithm spends its time. To this end we generated
larger problems, with the same structure, and we have broken down the timing
information. In Table 4 we consider a representative positive definite problem with "av-
erage characteristics", i.e., vary the problem
dimension n. (The sparsity structure remains the same.) The second column, labelled
"it", records the number of iterations required to achieve the stopping criteria; "totM"
records the total number of flops used by the (partial) Cholesky factorization ("m"
represents a million); "totls" records the number of flops used in the approximate line
search algorithm. Over 95% of the total flop count on these problems is represented by
the sum of the "totM" and "totls" columns - the remaining work in the algorithm, such
as the 2-dimensional trust region solution, is negligible in comparison (see Appendix
for more detail on the solution of trust region problems).
Observations on Table 4: First, there is no significant growth in number of iterations
as the problem dimension n increases. High accuracy is maintained for larger values
of n as well. As n increases the sparse matrix factorization work, "totM", increases
relative to the lines search cost, "totls". Therefore, speedup of the (partial) sparse
factorization aspect of the algorithm (e.g., use of parallelism, exploitation of
specific particular structure) will have significant impact on the overall computing time.
Conversely, improving the approximate line search (in terms of cost) is not a crucial
computing issue, at this point, for large-scale problems.
In addition to these randomly generated, but structured, positive definite problems,
we have experimented with three specific test cases. Two of these problems are from
the literature (e.g., [12, 24]) and the third example is new. In Tables 5 and 6 we report
on the "obstacle" problem - in the first case there are lower bounds only, in the second
case there are lower and upper bounds. In defining the specific example used we have
chosen the same parameter settings and specific functions used in [24]. Table 7 reports
on the elastic-plastic torsion problem. Again we used the same parameters as reported
in [24] to define the problem.
In

Table

8 we report on a linear spline approximation problem. This type of problem
arises, for example, in a particle method approach to turbulent combustion simulation
[28]. The problem results in a large sparse least-squares problem subject to
nonnegativity constraints on the variables. To set up a sample problem we assume
an m-by-m-by-m 3-dimensional grid. Within each cell are a set of particles randomly
located (we use approximately 10 particles per cell in our experiments). Each particle
p has a known function value, OE(p). Associate with each grid intersection point a linear
basis function and determine the best set of coefficients, x, for the basis functions, in
the least-squares sense, subject to nonnegativity constraints on x. The function OE we
used in our experiments is defined: given a point in 3-space,

Table
Obstacle Problem: Lower Bounds Only
its norm
50 2500 14 13
100 10,000 15 12

Table
Obstacle Problem: Lower and Upper Bounds
its norm
50 2500 14 12
100

Table
Elastic-plastic Torsion Problem
its norm
100 10,000 11 12
Observations on Tables 5-8. The most noteworthy observation is the apparent
insensitivity of our method to problem size for each of these problems. The number of
Linear Spline Approximation
its norm
22 10,648 17 11
iterations does not grow, for a given problem class, as the dimension of the problem
increases. For example, for the linear spline problem, 16 iterations are required when
iterations are required when Moreover, the number of
iterations is always modest, on this test set, i.e., less than 20. High accuracy is achieved
in all cases.
4.2. Indefinite Problems. We have adapted the Mor'e/Toraldo QP generation
scheme, in combination with sparse matrix functions in Matlab 4.0, to generate large
sparse indefinite matrices with a given sparsity pattern and given approximate set of
approximate eigenvalues. In the indefinite case we chose finite upper and lower bound
vectors, 1. (This is to avoid the generation of unbounded problems.)
In each of the problems in Tables 9-12 roughly 10% of the eigenvalues of H are neg-
ative. The column labels are the same as before however here "acc" does not represent
the number of accurate digits compared with the true solution since the true solution
is unknown due to indefiniteness of H. Instead, "acc" records the number of matching
digits in the objective function q in the last 2 iterations. (In each case the optimality
conditions were verified to hold at the final point.)

Table
Indefinite problems.
deg cond max avg acc
9 3 23 19.3 15
9 6 26 21.7 15
9 9
Observations on Tables 9-11: Iteration counts indicate that our method is not quite
as consistent or efficient on indefinite problems compared to the performance on positive
definite problems. Still, the overall efficiency seems very good - the average number of
Indefinite Problems.
deg cond max avg acc
9 3
9 6 19 17.7 15
3 9 14 11.3 15
9 9 25 18.3 15

Table
Indefinite problems.
deg cond max avg acc
9 3
9 6
6 9 15 13.
iterations required for any problem category is always less than 23.
In

Table

12 we indicate where the algorithm spends its time on indefinite problems
by considering a representative example and increasing the dimension.

Table
Indefinite Problems: Timing Breakdown
Remark on Table 12: We see no apparent growth in required iterations as n increases.
Clearly the "totM" column dominates the "totls" column as n increases. Recall that
"totM" represents the matrix factorization flop count while "totls" represents the number
of total flops required by the line search procedure. Therefore, to obtain further
improvements in efficiency for this type of approach it is best to focus on the matrix
factorization aspect of the overall procedure.
5. Theory and Conclusions. The numerical results obtained to date strongly
support the notion that a reflective Newton method represents an efficient way to accurately
locate local minimizers of large-scale (indefinite) quadratic functions subject to
bounds on some of the variables. The theory is supporting also: our reflective Newton
method is globally and quadratically convergent. Coleman and Li [10] present important
theoretical properties of reflective Newton methods for general nonlinear functions,
subject to bounds on some of the variables. The method in this paper is a specialization
of the general method to the quadratic case. Therefore, the general theory applies.
We make a compactness assumption before formally stating the main result.
Compactness Assumption: Given initial point x 1 2 F , it is assumed that the level
set
Theorem 2. Let fx k g be generated by Algorithm 5 with fs k g generated by Algorithm
6 and with fff k g determined by the approximate line search algorithm (Algorithm
7). If -
ffl Every limit point of fx k g is a first-order point.
ae M is the maximum spectral radius of -
M (x) on
)g. Since ae( -
is continuous on L, a compact set, the upper bound ae M exists. Recall that - 3
is a constant used in
Algorithm 6.
ffl Every nondegenerate limit point satisfies the second-order necessary conditions.
ffl If a nondegenerate limit point x   satisfies second-order sufficiency conditions
is sufficiently small, fx k g is convergent to x   ; the convergence
rate is quadratic, i.e.,
Proof. This algorithm is in the class of algorithm described in [10] and all the
assumptions of Theorem 20 in [10] are satisfied. The result follows from Theorem 20 in
[10].
In conclusion, strong theoretical and computational results indicate that a reflective
Newton method is an efficient and reliable way to solve problem (1.1) to high accuracy.
The computational results reported in this paper support this claim.
6.

Acknowledgement

. We thank our colleagues Jianguo Liu and Danny Ralph
for many helpful discussions on this work. Danny Ralph drew our attention to the
topology reference [25].
7.

Appendix

: The trust region problem. The trust region problem is
where A is a real symmetric matrix and k \Delta k denotes the 2-norm. The purpose of this
section is to review the nature of problem (7.1) and discuss a possible solution suitable
for low-dimensional problems. (In the context of our reflective Newton method for
problem (1.1), A is matrix -
M(x), a symmetric matix of order 2. The computational cost
of the procedure we describe to solve (7.1) is negligible compared to the other required
computations in the reflective Newton algorithm we propose.) For larger problems a
more approximate procedure is usually preferred, e.g., [14, 23, 30, 31]. Much of the
material in this section can be found elsewhere, e.g., [2, 4, 11, 14, 23, 29, 30, 31].
Diagonalization. We begin with an extremely useful characterization of the global
solution to (7.1).
Theorem 3. Vector s solves (7.1) if and only if there exists a scalar - 0 such
that
(a)
positive semidefinite;
(c) ksk - \Delta;
Proof. A proof is given, for example, in [30].
The usefulness of this result is best revealed after diagonalization. Suppose
where the columns of V are the orthonormal eigenvectors of A and
Obviously then so (a) is equivalent to
s:
By (b), - \Gamma- 1 , and so all vectors s satisfying (a, b) are of the form
The vector fi is arbitrary with respect to (a,b) but can sometimes help with respect
to satisfying (c,d). A basis for an algorithmic approach to this problem is to assume
the form given in (7.4) and strive to satisfy (c,d) by choosing -, and in some cases
fi, appropriately. (fi plays a role only if -  is the value of - at the
solution.)
The situation where (c,d) can be satisfied with easily
dispensed with (first half of Case 1 below). Therefore the primary task, assuming form
(7.4), is to determine -, and sometimes fi, to satisfy
We divide our approach into three possibilities. g.
Case 1: In this case either the Newton step is within the sphere ksk 2 - \Delta or it
is not. If kA then the optimal solution s
Case 2:

Figure

Obviously ks(-)k !1 as - ! \Gamma-
1 , and
1. Moreover, ks(-)k is convex; therefore, ks(-)k intersects \Delta in
exactly one place for - ? \Gamma- 1 .
Case 3:
ks 1 (\Gamma- 1 )k
is finite. Consider figures 3,4. There are now two possibilities. If ks 1 (\Gamma- 1
then there is a solution to (7.5) to the right of \Gamma- 1 ,
chosen to ensure (7.5), and -  . Note that if jIj ? 1
then the null space component of s may not be unique.
The Reciprocal Secular Equation. In theory we can build an algorithm on the
remarks given above. However it is better, numerically, to replace condition (7.5) with
\Gammaks(-)k
0:
Equation (7.6) is more linear in shape than equation (7.5); therefore, equation (7.6) is
more amenable to solution via Newton's method.
Considering the definition of rsec(-),
-n+-
it is easy to verify the following:
1. rsec(-) is convex on (\Gamma- 1 ; 1] and lim -!1
2. lim -!\Gamma- +rsec(-) is finite.
3. If ff i 6= 0, for some i 2 I, then rsec(\Gamma-
. Obviously
a single zero of rsec exits to the right of \Gamma- 1 in this case.
4. If 8i 2 I; ff
Algorithmic and Numerical Concerns. We assume that a solution to (7.1) is
sought and we are willing and able to compute full eigenvalue information,
(If this is not the case, perhaps due to the cost, then it is possible to approximately
solve (7.6) using an iterative scheme involving the Cholesky factorization of A
[14,
The method using appears to be straightforward.
Case 1: then the Newton step, \GammaA \Gamma1 g, is the solution with
then we can determine the zero of rsec(- ? 0.
Case 2:
and rsec admits
a single solution to the right of \Gamma- 1 .
Case 3:
If ks 1 (\Gamma- 1 )k - \Delta then there is a zero of rsec(-) to the right of \Gamma- 1 with
Otherwise, a solution to (7.1) is given by (7.4),
ks
Unfortunately, the situation is not quite so clean from a numerical point of view:
there is fuzziness between the second and third cases. In particular, if ff 1 is small the
equation very ill-conditioned for - near \Gamma- 1 and it can be quite difficult
(nigh impossible) to compute - such that rsec(-) is small. This extreme ill-conditioning
is due to the following "disagreement".
Assume I = f1g, for simplicity, and note that, if ff
On the other hand, if
which is not, in
general, the limiting value of (7.8). Therefore nearby problems (ff
can yield very different solutions to (7.6), and this is the cause of the ill-conditioning
of (7.6). Our solution to this ill-conditioning problem (trust.m) is to first attempt to
find a solution to However, if jrsec( -)j is not small, where -
- is
the computed "zero" returned by the zero-finder, then we set -
solution to (7.1) via (7.4).
This strategy works because the solution to (7.6) with ff i small for i 2 I is close to
a solution of (7.1) with the corresponding ff i at zero. To see this, initially assume that
Now first consider the case where ff the solution to (7.1) is:
where
otherwise, the result
is obvious). Define s
the solution to (7.1) can be written
Next consider ff small number. We can write the solution to (7.1) as
But as
which implies
Therefore, the solution to (7.1) with ff is near to the solution to (7.1) with ff
In general, if jIj ? 1 a solution to (7.1), with some components ff i near to zero, is
near to a solution of (7.1) with those components set to zero. In this case, where several
coefficients ff i equal to zero, i 2 I, problem (7.1) does not enjoy a unique solution [see
(7.4)] but the range space component is unique.



--R


Approximate solution of the trust region problem by minimization over two-dimensional subspaces
Advanced Computing Research Institute
Computing a trust region step for a penalty function
A direct active set algorithm for large sparse quadratic programs with simple bounds

A quadratically-convergent algorithm for the linear programming problem with lower and upper bounds



Global convergence of a class of trust region algorithms for optimization with simple bounds
On the minimization of quadratic functions subject to box constraints
Minimization of a quadratic function of many variables subject only to lower and upper bounds
Computing optimal locally constrained steps
Sparse matrices in matlab: Design and implemen- tation
Minimization subject to bounds on the variables

Curvilinear path steplength algorithms for minimization algorithms which use directions of negative curvature

A globally convergent method for l p problems
Solving the minimal least squares problem subject to bounds on the variables




A generalized conjugate gradient algorithm for solving a class of quadratic programming problems

Application of the velocity-dissipation PDF model to inhomogeneous turbulent flows
A family of trust-region-based algorithms for unconstrained minimization with strong global convergence properties
Trust region methods for unconstrained optimization
The conjugate gradient methods and trust regions in large scale optimization
A class of methods for solving large convex quadratic programs subject to box constraints
--TR

--CTR
Kamin Whitehouse , David Culler, Calibration as parameter estimation in sensor networks, Proceedings of the 1st ACM international workshop on Wireless sensor networks and applications, September 28-28, 2002, Atlanta, Georgia, USA
Kamin Whitehouse , David Culler, Macro-calibration in sensor/actuator networks, Mobile Networks and Applications, v.8 n.4, p.463-472, August
D. C. Jamrog , R. A. Tapia , Y. Zhang, Comparison of two sets of first-order conditions as bases of interior-point Newton methods for optimization with simple bounds, Journal of Optimization Theory and Applications, v.113 n.1, p.21-40, April 2002
Keiji Yanai , Nikhil V. Shirahatti , Prasad Gabbur , Kobus Barnard, Evaluation strategies for image understanding and retrieval, Proceedings of the 7th ACM SIGMM international workshop on Multimedia information retrieval, November 10-11, 2005, Hilton, Singapore
V. G. Domrachev , O. M. Poleshuk, A Regression Model for Fuzzy Initial Data, Automation and Remote Control, v.64 n.11, p.1715-1723, November
interior-point method;interior Newton method;quadratic programming