Tensor Methods for Large, Sparse Unconstrained Optimization.
Tensor methods for unconstrained optimization were first introduced by Schnabel and Chow [SIAM J. Optim., 1 (1991), pp. 293--315], who described these methods for small- to moderate-sized problems. The major contribution of this paper is the extension of these methods to large, sparse unconstrained optimization problems. This extension requires an entirely new way of solving the tensor model that makes the methods suitable for solving large, sparse optimization problems efficiently. We present test results for sets of problems where the Hessian at the minimizer is nonsingular and where it is singular. These results show that tensor methods are significantly more efficient and more reliable than standard methods based on Newton's method.
Introduction
In this paper we describe tensor methods for solving the unconstrained optimization problem
where D is some open set containing x   . We assume that f is at least twice continuously
differentiable, and n is large.
Tensor methods for unconstrained optimization are general purpose methods primarily intended
to improve upon the performance of standard methods, especially on problems where
deficiency. They are also intended to be at least as efficient as standard
methods on problems where r 2 f(x   ) is nonsingular.
Tensor methods for unconstrained optimization base each iteration upon the fourth order
model of the objective function f(x)
is the current iterate, rf(x c ) and r 2 f(x c ) are the first and second analytic
derivatives of f at x c , or finite difference approximations to them, and where the tensor terms
at x c , T c 2 ! n\Thetan\Thetan and V c 2 ! n\Thetan\Thetan\Thetan , are symmetric. (We use the notation rf(x c ) \Delta d for
)d to be consistent with the tensor notation T c \Delta d 3 and
Also, for simplicity, we abbreviate terms of the form dd; ddd, and dddd by d 2 ; d 3 , and d 4 ,
respectively.) Before proceeding, we define the tensor notation used above.
n\Thetan\Thetan . Then for
n\Thetan\Thetan\Thetan . Then for
The tensor terms are selected so that the model interpolates a small number of function and
gradient values from previous iterations. This results in T c and V c being low-rank tensors, which
is crucial for the efficiency of the tensor method. The tensor method requires no more function
or derivative evaluations per iteration and hardly more storage or arithmetic operations, than a
standard method based on Newton's method.
Standard methods for solving unconstrained optimization problems are widely described in
the literature; general references on this topic include Dennis and Schnabel [9], Fletcher [11],
and Gill, Murray, and Wright [13]. In this paper, we propose extensions to standard methods
that use analytic or finite difference gradients and Hessians.
The standard method for unconstrained optimization, Newton's method, bases each iteration
upon the quadratic model of f(x)
This method is defined when r 2 f(x c ) is nonsingular, and consists of setting the next iterate x+
to the minimizer of (1.3), i.e.,
A distinguishing feature of Newton's method is that if r 2 f(x c ) is nonsingular at a local
minimizer x   , then the sequence of iterates produced by (1.4) converges quadratically to x   .
However, Newton's method is generally linearly convergent at best if r 2 f(x   ) is singular [14].
Methods based on (1.2) have been shown to be more reliable and more efficient than standard
methods on small to moderate size problems [18]. In the test results obtained for both non-singular
and singular problems, the improvement by the tensor method over Newton's method
is substantial, ranging from 30% to 50% in iterations, and function and derivative evaluations.
The improvement is even more dramatic for singular problems. Furthermore, the tensor method
solves several problems that Newton's method fails to solve.
The tensor algorithms described in [18] are QR-based algorithms involving orthogonal transformations
of the variable space. These algorithms are very effective for minimizing the tensor
model when the Hessian is dense because they are very stable numerically, especially when the
Hessian is singular. However, they are not efficient for sparse problems because they destroy the
sparsity of the Hessian due to the orthogonal transformation of the variable space. To preserve
the sparsity of the Hessian, we have developed an entirely new way of solving the tensor model
that employs a sparse variant of the Cholesky decomposition. This makes our new algorithms
very well suited for sparse problems.
The remainder of this paper is organized as follows. In x2 we briefly review the techniques
used to form the tensor model, that were introduced in Schnabel and Chow [18]. In x3 we describe
efficient algorithms for minimizing the tensor model when the Hessian is sparse. xx4 and
5 discuss the globally convergent modifications for tensor methods for unconstrained optimiza-
tion. These consist of line search backtracking and model trust region techniques. A high level
implementation of the tensor method is given in x6. In x7 we describe comparative testing for
an implementation based on the tensor method versus an implementation based on Newton's
method, and present summary statistics of the test results. Finally, a summary of our work and
a discussion of future research is given in x8.
2. Forming the Tensor Model
In this section, we briefly review the techniques for forming the tensor model for unconstrained
optimization that were introduced in [18].
As was stated in the previous section, the tensor method for unconstrained optimization
bases each iteration upon the fourth order model of the nonlinear function f(x) given by (1.2).
The choices of T c and V c in (1.2) cause the third order term T c \Delta d 3 and the fourth order
to have simple and useful forms. These tensor terms are selected so that the tensor
model interpolates function and gradient information at a set of p not necessarily consecutive
past iterates x
In the remainder of this paper, we restrict our attention to 1. The reasons for this
choice are that the performance of the tensor version that allows p - 1 is similar overall to that
constraining p to be 1, and that the method is simpler and less expensive to implement in this
case. (The derivation of the third and fourth order tensor terms for p - 1 is explained in detail
in [18].)
The interpolation conditions at the past point x are given by
and
where
Schnabel and Chow [18] choose T c and V c to satisfy (2.1) and (2.2). They first show that
the interpolation conditions (2.1) and (2.2) uniquely determine T c \Delta s 3 and V c \Delta s 4 . Multiplying
(2.2) by s yields
Let ff, fi 2 ! be defined by
Then from (2.1) and (2.3) they obtain the following system of two linear equations in the two
unknowns ff and fi:2 ff
are defined by
The system (2.4)-(2.5) is nonsingular, therefore the values of ff and fi are uniquely determined.
Hence, the interpolation conditions uniquely determine T c \Delta s 3 and V c \Delta s 4 . Since these are the
only interpolation conditions, the choice of T c and V c is vastly underdetermined.
Schnabel and Chow [18] choose T c and V c by first selecting the smallest symmetric V c , in the
Frobenius norm, for which
where fi is determined by (2.4)-(2.5). Then they substitute this value of V c into (2.2), obtaining
where
This is a set of n linear equations in n 3 unknowns T c (i; j; k), 1 - Finally, Schnabel
and Chow [18] choose the smallest symmetric T c and V c , in the Frobenius norm, which satisfy
the equations (2.6)-(2.7). That is,
min
Vc2! n\Thetan\Thetan\Thetan
subject to V c \Delta s
and
min
n\Thetan\Thetan
subject to T c \Delta s
The solution to (2.8) is
(s\Omega s\Omega s\Omega s);
where the tensor V
n\Thetan\Thetan\Thetan is called a fourth order rank-one tensor for
which use the
notation\Omega to be consistent
with [18].)
The solution to (2.9) is
s\Omega s\Omega b; (2:10)
where the notation
n\Thetan\Thetan , is called a third order rank-one
tensor for which T (i; j; is the unique vector for which (2.10) satisfies
(2.6), and is given by
determined by the minimum norm problems (2.9) and (2.8) have rank 2 and 1,
respectively. This is the key to form, store, and solve the tensor model efficiently.
The whole process of forming the tensor model requires only O(n 2 ) arithmetic operations. The
storage needed for forming and storing the tensor model is only a total of 6n.
For further information we refer to [18].
3. Solving the Tensor Model when the Hessian is Sparse
In this section we give efficient algorithms for finding a minimizer of the tensor model (1.2),
when the Hessian is sparse.
The substitution of the values of T c and V c into (1.2) results in the tensor model
As we stated in x2, we only consider the case where the tensor model interpolates f(x) and
rf(x) at the previous iterate, i.e., 1. The generalization for p - 1 is fairly straightforward.
This constraint is mainly motivated by our computational results. When we allow p - 1, our
test results showed almost no improvement over the case where 1. The tensor method is
therefore considerably simpler, and cheaper in terms of storage and cost per iteration.
3.1. Case 1: the Hessian is Nonsingular
We show that the minimization of (3.1) can be reduced to the solution of a third order polynomial
in one unknown, plus the solution of three systems of linear equations that all involve the same
coefficient matrix r 2 f(x c ). For conciseness, we use the notation
A necessary condition for d to be a local minimizer of (3.1) is that the derivative of the
tensor model with respect to d must be zero. That is,
which yields
If we first premultiply equation (3.2) by s T on both sides, we obtain a cubic equation (in fi) in
the unknowns
If we then premultiply equation (3.2) by b T on both sides, we obtain another cubic equation (in
fi) in the unknowns fi and ',
Thus, we obtain a system of two cubic equations in the two unknowns fi and ' which can be
solved analytically.
We now show how to compute the solutions of this system of two cubic equations in two
unknowns by computing the solutions of a single cubic equation in the unknown fi. Let
first calculate the
value of ' as a function of fi using equation (3.3), i.e.,
Note that the denominator of equation (3.5) is equal to zero if either
assume that fi 6= 0, otherwise the tensor model would be reduced to the Newton model. Now,
would be quadratic in fi, therefore
Thus, real valued minimizers of the tensor model (3.1) may exist only if 0: It is easy
to check that in order for ' to have a defined cannot be zero.
If fi 6= 0 and w 6= 0, we substitute expression for ' into equation (3.4) and obtain
which is a third order polynomial in the one unknown fi. The roots of equation (3.6) are
computed analytically. We substitute the values of fi into equation (3.5) to calculate the values
of '. Then we simply substitute the values of fi and ' into equation (3.2) to obtain the values
of d. The major cost in this whole process is the calculation of H
After we compute the values of d, we determine which of them are potential minimizers.
Our criterion is to select those values of d which guarantee that there is a descent path from x c
to x c + d for the model M T among the selected steps, we choose the one that is
closest to the current iterate x c in the euclidean norm sens. If the tensor model has no minimizer
we use the standard Newton step as the step direction for the current iteration.
3.2. Case 2: the Hessian is Rank Deficient
If the Hessian matrix is rank deficient we transform the tensor model given in (3.1) by the
following procedure. Let
d, and ffi is the new unknown. Substituting this
expression for d into (3.1) yields the following tensor model which is a function of ffi,
d) 2
d)
d)s
d)
(b T -
If we let -
d, -
d, -
we obtain the modified tensor model,
d)
The advantage of this transformation is that the matrix -
H is likely to be nonsingular if the rank
of (r 2 f(x c )) is at least n \Gamma 1. A necessary and sufficient condition for -
H to be nonsingular is
given in the following lemma. Let g and H denote rf(x c ) and r 2 f(x c ), respectively.
Lemma 3.1. Let H 2 ! n\Thetan , s
css T is nonsingular if and only if M =6 6 6 6 6 4
H cs
cs
is nonsingular:
(Note that the
submatrix was premultiplied by the constant c to symmetrize the
augmented matrix M .)
Proof. We prove that there exists only if there
exist
H cs
cs
Suppose first that (H
Conversely, if there exists (-v; w) satisfying (3.9), then s T -
otherwise, contradicts (3.9). Thus (H singular if and only if M is
singular.
Corollary 3.2. Let H 2 ! n\Thetan , s
css T is nonsingular then
H cs
has full row rank
Proof. follows from Lemma 3.1.
Lemma 3.3. Let H 2 ! n\Thetan ,
css T is nonsingular if and only if
H cs
has full row rank.
Proof. The only if part follows from Corollary 3.2. Now assume
H cs
has full row rank.
Since H has rank
n\Theta(n\Gamma1) have full column rank. Since
H cs
has full row rank,
(v
From
T and H 2 has full column rank, (3.10) is equivalent to
(v
Thus the n \Theta n matrix
is nonsingular. Analogously, the n \Theta n matrix
is
nonsingular. Therefore
is nonsingular. 2
For ffi to be a local minimizer of (3.8) the derivative of the tensor model (3.8) with respect
to ffi must be zero. That is,
which yields
Premultiplying equation (3.12) by s T on both sides results in a cubic equation (in fi) in the two
unknowns
fis T -
fis T -
fis T -
The premultiplication of equation (3.12) by b T on both sides yields another cubic equation (in
fi) in the two unknowns fi and '
Therefore, we obtain a system of two cubic equations in the two unknowns fi and ' which we
can solve analytically.
Since equation (3.13) is linear in ', we can compute ' as a function of fi, and then substitute
its expression into equation (3.14) to obtain an equation in the one unknown fi. Let
g, and
The denominator of equation (3.15) is equal to zero if either -
then equation (3.13) would be quadratic in fi, therefore
Hence, real valued minimizers of the tensor model (3.8) may exist only if (1
It is straightforward to verify from (3.14) that for ' to be defined
reduces to the following cubic equation in fi
Once we calculated the expressions for fi from equation (3.16), we substitute them into the
following equation for ' obtained from equation (3.14)
If neither -
equation (3.14) and obtain
fiw \Gamma2
which is a third order polynomial in the one unknown fi. The roots of equation (3.17) are then
computed analytically. After we determine the values of fi, we substitute them into equation
(3.15) to calculate the corresponding values of '. then, we simply substitute the values of fi
and ' into equation (3.12) to obtain the values of ffi . The dominant cost in this whole process is
the computation of -
Similar to the nonsingular case, a minimizer ffi is selected such that there exists a descent
path from the current point x c to x c that it is closest to x c .
To obtain the tensor step d we set d to -
An appropriate choice of -
d is the step used in
the previous iteration simply because it has the right scale.
To solve linear systems of the form -
n\Thetan sparse
and we use the augmented matrix M defined in Lemma 3.1. That is, we write
H cs
cs
x
The (n in (3.18) is sparse and can be factorized efficiently as long as the
last row and column are not pivoted until the last few iterations. In fact, we can combine the
nonsingular and singular cases by factorizing H , but we shift to a factorization of the augmented
matrix if H is discovered to be singular with rank n \Gamma 1. However, we use a Schur complement
method to obtain the solution of the augmented matrix by updating the solution from the system
b. This choice was motivated by the fact that the Schur complement method was simpler
and more convenient to use than the factorization of the augmented matrix M . We describe
this updating scheme in x6.
If the Schur complement method shows that M is rank deficient (a case that is very rare in
practice), or H has rank less than use the standard Newton step as the step direction
for the current iteration.
4. Line Search Backtracking Techniques
The line search global strategy we used in conjunction with our tensor method for large sparse
unconstrained optimization is similar to the one used for nonlinear equations [4, 6]. This strategy
has shown to be very successful for large sparse systems of nonlinear equations. We also
found that it is superior to the approach used by Schnabel and Chow [18]. The main difference
between the two approaches is that ours always tries the full tensor step first. If this provides
enough decrease in the objective function then we terminate, otherwise we find acceptable next
iterates in both the Newton and tensor directions and select the one with the lower function
value as the next iterate. Schnabel and Chow on the other hand, always find acceptable next
iterates in both the Newton and tensor directions and choose the one with the lower function
value as the next iterate. In practice, our approach almost always requires fewer function evaluations
while retaining the same efficiency in iteration numbers. The global framework for line
search methods for unconstrained minimization is given in Algorithm 4.1.
Algorithm 4.1. Global Framework for Line Search Methods for Unconstrained Minimization.
Let x c be the current iterate
d t the tensor step
d n is the Newton step
and
if (minimizer of the tensor model was found) then
else
Find an acceptable x n
in the Newton direction d n
using Algorithm A6.3.1 page 325 ([9])
Find an acceptable x t
in the tensor direction d t
using Algorithm A6.3.1 page 325 ([9])
if f(x n
else
endif
endif
else
Find an acceptable x n
in the Newton direction d n
using Algorithm A6.3.1 page 325 ([9])
endif
5. Model Trust Region Techniques
The two computational methods that are generally used for approximately solving the trust
region problem based on the standard model,
subject to jj d jj
where ffi c is the current trust region radius, are the locally constrained optimal (or "hook") step,
and the dogleg step. When ffi c is shorter than the Newton step, the locally constrained optimal
step [16] finds a - c such that jj d(- c
takes The dogleg step is a modification of the trust region algorithm introduced
by Powell [17]. However, rather than finding a point on the curve d(- c ) such
that it approximates this curve by a piecewise linear function in the subspace
spanned by the Newton step and the steepest descent direction \Gammarf (x c ), and takes x+ as the
point on this approximation for which jj x+ \Gamma x c e.g. [9] for more details.)
Unfortunately these two methods are hard to extend to the tensor model, which is a fourth
order model. Trust region algorithms based on (5.19) are well defined because it is always
possible to find a unique point x+ on the curve such that jj x+ \Gamma x c . Additionally, the
value of f(x c )+rf(x c ) along the curve d(- c ) is monotonically decreasing from
x c to x n
which makes the process reasonable. These properties do not
extend to the tensor model which is a fourth order model that may not be convex. Furthermore,
the analogous curve to d(- c ) is more expensive to compute. For these reasons, we consider a
different trust region approach for our tensor methods.
The trust region approach that is discussed in this section is a two-dimensional trust region
step over the subspace spanned by the steepest descent direction and the tensor (or standard)
step. The main reasons that lead us to adopt this approach is because it is easy to construct,
closely related to dogleg type algorithms over the same subspace. This step may be close to
optimal trust region step algorithms in practice. Byrd, Schnabel, and Shultz [7] have shown that
for unconstrained optimization using a standard quadratic model, the analogous two-dimensional
minimization approach produces nearly as much decrease in the quadratic model as the optimal
trust region step in almost all cases.
The two-dimensional trust region approach for the tensor model computes an approximate
solution to
subject to jj d jj
by performing a two-dimensional minimization,
subject to jj d jj
where d t and g s are the tensor step and the steepest descent direction, respectively, and ffi c is the
trust region radius. This approach will always produce a step that reduces the quadratic model
by at least as much as a dogleg type algorithm which reduces d to a piecewise linear curve in the
same subspace. At each iteration of the tensor algorithm, the trust region method either solves
(5.20), or minimizes the standard linear model over the two-dimensional subspace spanned by
the standard Newton step and the steepest descent direction. The decision of whether to use
the tensor or standard model is made using the following criterion:
if (no minimizer of the tensor model was found) or (rf(x c
then
selected by trust region algorithm
else
selected by trust region algorithm
endif
Before we define the two-dimensional trust region step for tensor methods, we show how to
convert the problem
subject to jj d jj
to an unconstrained minimization problem.
First, we make g s orthogonal to d t by performing the Householder transformation:
then, we normalize both - g s and d t to obtain:
~
~
s
Since d is in the subspace spanned by the tensor step ~
d t and the steepest descent direction ~
s ,
it can be written as
If we square the l 2 norm of this expression for d and set it to ffi 2 , we obtain the following equation
for fi as a function of ff
Substituting this expression for fi into (5.25) and then the resulting d into (5.21), yields the
global minimization problem in the one variable ff, given by (5.26) bellow. Thus, problems
and (5.21) are equivalent. Let g hg = ~
s .
c
)ff
To transform the problem
subject to jj d jj
to an unconstrained minimization problem, we use the same procedure described above to show
that (5.27) is equivalent to the following global minimization problem in the one variable ff:
c
Algorithm 5.1. Two-Dimensional Trust Region for Tensor Methods
Let d t be the tensor step
d n the standard step
x c the current iterate
x+ the next iterate
steepest descent direction
and ffi c the current trust region radius.
~
are given by (5.23) and (5.24), respectively.
~
d n is obtained in an analogous way to ~
applying transformations (5.22) and (5.23) to it.
1. if tensor model selected then
Solve problem (5.26) using the procedure described in Algorithm 3.4 [6]
else fstandard Newton model selectedg
Solve problem (5.28) using the procedure described in Algorithm 3.4 [6]
endif
2. if tensor model selected then
~
s

where ff   is the global minimizer of (5.26)
else fstandard Newton model selectedg
~
s

where ff   is the global minimizer of (5.28)
endif
3. f Check new iterate and update trust region radius.g
pred
the global step d is successful
else
decrease trust region
go to step 1
endif
where
pred
pred
standard Newton model selected.
The methods used for adjusting the trust radius during and between steps, are given in Algorithm
page 338 ([9]). The initial trust radius can be supplied by the user, if not, it is set to the
length of the initial Cauchy step.
6. A High Level Algorithm for the Tensor Method
In this section, we present the overall algorithm for the tensor method for large sparse unconstrained
optimization. Algorithm 6.1 is a high level description of an iteration of the tensor
method, that was described in xx3-5. A summary of the test results for this implementation is
presented in x7.
Algorithm 6.1. An Iteration of the Tensor Method for Large Sparse Unconstrained Optimization

Let x c be the current iterate
d t the tensor step
and d n the Newton step.
1. Calculate rf(x c ) and decide whether to stop. If not:
2. Calculate r 2 f(x c ).
3. Calculate the terms T c and V c in the tensor model, so that the tensor model interpolates
f(x) and rf(x) at the past point.
4. Find a potential minimizer d t of the tensor model (3.1).
5. Find an acceptable next iterate x+ using either a line search or a two-dimensional trust
region global strategy.
go to step 1.
In step 1, the gradient is either computed analytically or approximated by the algorithm
A5.6.3 given in Dennis and Schnabel [9]. In step 2, the Hessian matrix is either calculated
analytically or approximated by a graph coloring algorithm described in [8]. Note that it is
crucial to supply an analytic gradient if the finite difference Hessian matrix requires many
gradient evaluations. Otherwise, the methods described in this paper may not be practical, and
inexact type of methods may be preferable. The procedures for calculating T c and V c in step 3
were discussed in x2. Step 4 calculates d t as described in xx 3-4. The Newton step d n is also
computed as a by product of the minimization of the tensor model. The Newton step d n is
the modified Newton step (r 2 f(x c safely positive
definite, and - ? 0 otherwise. To obtain the perturbation - we use a modification of MA27 [10]
advocated by Gill, Murray, Ponceleon, and Saunders in [12]. In this method we first compute the
LDL T of the Hessian matrix using the MA27 package, then change the block diagonal matrix
D to D + E. The modified matrix is block diagonal positive definite. This guarantees that the
E)L T is positive definite as well. Note, that the Hessian matrix is not
modified if it is already positive definite.
The tensor and Newton algorithms terminate if jj rf(x c ) jj
Another implementation issue that deserves some attention is how to find a solution to the
augmented system (3.18), when the Hessian matrix is rank deficient. To do this, we use a Schur
complement method to update the solution x obtained from solving Hx = b. This requires
that H must have full rank. Thus, some modifications are necessary in order for this method
to work. We have modified the factorization phase of MA27 to be able to detect the row and
column indices of the first pivot that is less or equal than some given tolerance tol. Note that if
the rank of the Hessian matrix is less than we skip this whole updating scheme and
perturb the matrix as described in the previous paragraph. We also modified the solve phase
of MA27 such that whenever there is a zero pivot, the corresponding solution component is set
to zero. This way the solution of is the same as the solution of H e
is the matrix H minus the row and column at which singularity occurred. Since y has
components, the remaining one, which is also the component corresponding to the zero pivot, is
set to 0.) Afterwards, we obtain the solution of an augmented system using a Schur complement
method, where the coefficient matrix is the matrix H augmented by two rows and columns, i.e.,
the (n+ 1)-st row and column are the ones at which singularity was detected, and the (n+2)-nd
row and column are cs T and cs, respectively. The Schur complement method is implemented by
first invoking MA39AD [1] to form the Schur complement H in the extended
matrix, where D is the 2 by 2 lower right submatrix, C is the lower left 2 by n submatrix, and
B is the upper right n by 2 submatrix, of the augmented matrix. The Schur complement is then
factorized into its QR factors. Next, MA39BD [1] solves the extended system (3.18) using the
following well-known scheme:
1. Solve
2. Solve y.
3. Solve
4.
7. Test Results
We tested our tensor and Newton algorithms on a variety of nonsingular and singular test
problems. In the following we present and discuss summary statistics of the test results.
All our computations were performed on a SUN using double
precision arithmetic.
First, we tested our program on the set of unconstrained optimization problems from the
[3] and the MINPACK-2 [2] collections. Most of these problems have nonsingular Hessians
at the solution. We also created singular test problems as proposed in [4, 19] by modifying the
nonsingular test problems from the CUTE collection as follows. Let
be the function to minimize, where f is the number of element functions,
and
In many cases, F at the minimizer x   , and F 0 (x   ) is nonsingular. then according to
[4, 19], we can create singular systems of nonlinear equations from (7.1) by forming
n\Thetak has full column rank with 1 - k - n. Hence, -
k. For unconstrained optimization, we simply need to define the singular function
From (7.3) and -
and
we know that r 2 -
By using (7.2) and (7.3), we created two sets of singular problems, with r 2 -
respectively, by using
and
respectively. The reason for choosing unit vectors as columns for the matrix A is mainly to
preserve the sparsity of the Hessian during the transformation (7.2).
For all our test problems we used a standard line search backtracking strategy. All the test
problems with the exception of rank problems were ran with analytic
gradients and Hessians provided by the CUTE and MINPACK-2 collections. For rank
test problems, we have modified the analytic gradients provided by the CUTE collection
to take into account the modification (7.2). On the other hand, we used the graph coloring
algorithm [8] to evaluate the finite difference approximation of the Hessian matrix.
A summary for the test problems whose Hessians at the solution have ranks n,
presented in Table 1. The descriptions of the test problems and the detailed results are
given in the Appendix. In Table 1 columns "better" and "worse" represent the number of times
the tensor method was better and worse, respectively, than Newton's method by more than one
gradient evaluation. The "tie" column represents the number of times the tensor and standard
methods required within one gradient evaluation of each other. For each set of problems, we
summarize the comparative costs of the tensor and standard methods using average ratios of
three measures: gradient evaluations, function evaluations, and execution times. The average
gradient evaluation ratio (geval) is the total number of gradients evaluations required by the
tensor method, divided by the total number of gradients evaluations required by the standard
method on these problems. The same measure is used for the average function evaluation
(feval) and execution time (time) ratios. These average ratios include only problems that were
successfully solved by both methods. We have excluded all cases where the tensor and standard
methods converged to a different minimizer. However, the statistics for the "better", "worse",
and "tie" columns include the cases where only one of the two methods converges, and exclude
the cases where both methods do not converge. We also excluded problems requiring a number
of gradient evaluations less or equal than 3 by both methods. Finally, columns "t/s" and "s/t"
show the number of problems solved by the tensor method but not by the standard method
and the number of problems solved by the standard method but not by the tensor method,
respectively.
The improvement by the tensor method over the standard method on problems with rank
averaging 48% in function evaluations, 52% in gradient evaluations, and 59% in
execution times. This is due in part to the rate of convergence of the tensor method being faster
than that of Newton's method, which is known to be only linearly convergent with constant3
. On problems with rank the improvement by the tensor method over the standard
method is also substantial, averaging 30% in function evaluations, 37% in gradient evaluations,
and 34% in execution times. In the test results obtained for the nonsingular problems, the tensor
method is 9% worse than the standard method in function evaluations, but 31% and 33% better
in gradient evaluations and in execution times, respectively. The main reason for the tensor
method requiring on the average more function evaluations than the standard method is because
on some problems, the full tensor step does not provide sufficient decrease in the objective
function, and therefore the tensor method has to perform a line search in both the Newton
and tensor directions, which causes the number of function evaluations required by the tensor
method to be inflated. As a result, we intend to investigate other possible global frameworks for
line search methods that could potentially reduce the number of functions evaluations for the
tensor method.
To obtain an experimental indication of the local convergence behavior of the tensor and
Newton methods on problems where rank(r 2 f(x   examined the sequence of ratios
produced by the Newton and tensor methods on such problems. These ratios for a typical
problem are given in Table 2. In almost all cases the standard method exhibits local linear
convergence with constant near 2, which is consistent with the theoretical analysis. The local
convergence rate of the tensor method is faster with a typical final ratio of around 0.01. Whether
this is a superlinear convergence remains to be determined. We have done similar experiments for
problems with rank(r 2 f(x   and the tensor method did not show a faster-than-linear
convergence rate, because it did not have enough information since
The tensor method solved a total of four nonsingular problems, five rank
and 7 rank that Newton's method failed to solve. The reverse never occurred.
This clearly indicates that the tensor method is most likely to be more robust than Newton's
method.
The overall results show that having some extra information about the function and gradient
in the past step direction is quite useful to achieve the advantages of tensor methods.
8. Summary and Future Research
In this paper we presented efficient algorithms for solving large sparse unconstrained optimization
using tensor methods. We described new methods for minimizing the tensor model, that are
efficient for problems where the Hessian matrix is large and sparse. Implementations using these
tensor methods have been shown to be considerably more efficient especially on problems where
Rank tensor/standard pbs solved Average Ratio-Tensor/standard

Table

1: Summary of the CUTE and MINPACK-2 test problems using line search
Iteration (k) Standard method Tensor method
9 0.600 0.126
14 0.969
22 0.896
26 0.667
28 0.666

Table

2: Speed of convergence on the BRYBND problem with rank(r 2 f(x   modified
by (7.2), started from x 0 . The ratios in second and third columns are defined by
(7.
the Hessian matrix has a small rank deficiency at the solution. Typical gains over standard
Newton methods range from 40% to 50% in function and gradient evaluations, and in computer
time. The size and consistency of the efficiency gains indicate that the tensor method may be
preferable to Newton's method for solving large sparse unconstrained optimization problems
where analytic gradients and/or Hessians are available. To firmly establish such a conclusion,
additional testing is required, including test problems of very large size.
On sparse problems where the function or the gradient is expensive to evaluate, the finite
difference approximation of the Hessian matrix by the graph coloring algorithm [8] may be
very costly. Hence, Quasi-Newton methods may be preferable to use in this case. These are
methods, which involve low-rank corrections to a current approximate Hessian matrix. We are
currently attempting to extend our tensor methods to Quasi-Newton methods for large sparse
unconstrained minimization problems.
We also considered solving large sparse structured unconstrained optimization problems using
tensor methods. In this variant, we explore the possibility of using exact third and fourth
order derivative information. The calculation of these derivatives is simplified using the concept
of partial separability, a structure which has already proven to be useful when building quadratic
models for large scale nonlinear problems [15]. However, the calculation of the minimizer of this
exact tensor model is more problematic because we need to solve a sparse system of nonlinear
equations. An obvious approach to solve these equations is to use a Newton-like method. Such
a method is characterized by the approximation of the Jacobian used in the Newton process. A
simple idea is to use a fixed Jacobian at each step. This has the advantage that the Jacobian
will have already been obtained in the current tensor iteration. However, potential slow convergence
of such a scheme may make the cost of a tensor iteration prohibitive. We are currently
investigating other possible approaches such as a modified Newton's method in which the approximated
Jacobian matrix will incorporate more useful information, or an iterative method
such as a nonlinear GMRES. This work, a cooperation with Nick Gould [5], will be reported in
the near future.
We are almost done with the implementation and testing of the two-dimensional trust region
global strategy described in x5. This work will be reported in a forthcoming paper.
We are also implementing the algorithms discussed in this paper in a software package. This
package uses one past point in the formation of the tensor terms, which makes the additional
cost and storage of the tensor method over the standard method very small. The package will
be available soon.

Acknowledgments

. We thank Professor Bobby Schnabel for his suggestions on how to
minimize the tensor model when the Hessian is rank deficient, Nick Gould for discussing a number
of implementation issues, Ta-Tung Chow for reviewing the first draft of the paper, and my
CERFACS colleage Jacko Koster for his numerous suggestions.



--R


The Minpack-2 test problem col- lection
CUTE: Constrained and Unconstrained Testing Environment.
Solving large sparse systems of nonlinear equations and nonlinear least squares problems using tensor methods on sequential and parallel computers.
Tensor methods for large-scale unconstrained opti- mization
TENSOLVE: a software package for solving systems of nonlinear equations and nonlinear least squares problems using tensor methods.
Approximation solution of the trust region problem by minimization over two-dimensional subspaces
Estimating sparse hessian matrices.
Numerical methods for unconstrained optimization and nonlinear equations.
A set of Fortran subroutines for solving sparse symmetric sets of linear equations.
Practical method of optimization
Preconditioners for indefinite systems arising in optimization and nonlinear least squares problems.
Practical Optimization.
Analysis of Newton's method at irregular singularities.
On the unconstained optimization of partially separable functions.
The Levenberg-Marquardt algorithm: implementation and theory
A new algorithm for unconstrained optimization.
Tensor methods for unconstrained optimization using second derivatives.
Tensor methods for nonlinear equations.
--TR
large-scale optimization;tensor methods;unconstrained optimization;sparse problems;singular problems