%\documentclass[glov2, final]{svjour2}
\documentclass[runningheads]{svjour2}
%
\smartqed  % flush right qed marks, e.g. at end of proof
%
%\usepackage{graphicx}
%\usepackage{biograph}        % to allow for author biography at the end
%
%\usepackage{mathptmx}      % use Times fonts if available on your TeX system
%
% insert here the call for the packages your document requires
%\usepackage{latexsym}
% etc.
%
\def\eqnok#1{(\ref{#1})}
\newcommand{\h}{\hat}
\newcommand{\hX}{\hat X}
\newcommand{\N}{{\cal N}}
\newcommand{\nn}{\nonumber} 
\newcommand{\inte}[1]{\bar{#1}}

\newtheorem{assumption}{Assumption}
\newcommand{\beq}{\begin{equation}}
\newcommand{\eeq}{\end{equation}}
\newcommand{\beqa}{\begin{eqnarray}}
\newcommand{\eeqa}{\end{eqnarray}}
\newcommand{\beqas}{\begin{eqnarray*}}
\newcommand{\eeqas}{\end{eqnarray*}}
\newcommand{\bi}{\begin{itemize}}
\newcommand{\ei}{\end{itemize}}
\newcommand{\gap}{\hspace*{2em}}
\newcommand{\sgap}{\hspace*{1em}}
\newcommand{\mgap}{\vspace{.1in}} 
\newcommand{\ov}{\overline}
\newcommand{\vgap}{\vspace*{.1in}}
\newcommand{\diag}{\rm diag}
\newcommand{\argmin}{{\rm argmin}}
\newcommand{\Argmin}{{\rm Argmin}}  
\newcommand{\Ker}{{\rm Ker}}
\newcommand{\R}{\rm Im}
\newcommand{\half}{{\textstyle{1\over 2}}}
%\newcommand{\proof}{\par{\it Proof}. \ignorespaces}
%\newcommand{\endproof}{{\ \hfill \hbox{\vrule width1.0ex height1.0ex}\parfillskip 0pt}\par}
%\newcommand{\endremark}{{\ \hfill\hbox{\vrule width1.0ex height1.0ex}\parfillskip 0pt}\par}
%\newcommand{\QED}{\ifhmode\unskip\nobreak\fi\ifmmode\ifinner\else\hskip5pt\fi\fi
%  \hbox{\hskip5pt\vrule width5pt height5pt depth1.5pt\hskip1pt}}
\newcommand{\tD}{{\tilde D}}
\newcommand{\tE}{{\tilde E}}
\newcommand{\tV}{{\tilde V}}
\newcommand{\tp}{{\tilde p}}        
\newcommand{\Mh}{{M^{1/2}}}
\newcommand{\Mhn}{{M^{-1/2}}}
\newcommand{\Dh}{{D^{1/2}}}
\newcommand{\Dhn}{{D^{-1/2}}}
\newcommand{\lambdaL}{{\lambda^L}}
\newcommand{\lambdaU}{{\lambda^U}}
\newcommand{\tDn}{{{\tilde D}^{-1}}}
\newcommand{\tEn}{{{\tilde E}^{-1}}}
\newcommand{\tW}{{\tilde W}}
\newcommand{\hDn}{{{\hat D}^{-1}}}
\newcommand{\hEn}{{{\hat E}^{-1}}}
\newcommand{\hV}{{\hat V}}
\newcommand{\hW}{{\hat W}}  
\newcommand{\lambdan}{{\lambda_n}}
\newcommand{\bv}{{\bar V}}  
\newcommand{\bD}{{\bar D}}  
\newcommand{\rank}{\mbox{\rm rank}}
\newcommand{\Diag}{\mbox{\rm Diag}}
\newcommand{\cM}{{\cal M}}   
\newcommand{\lambdak}{{\lambda^{(k)}}}
\newcommand{\muk}{{\mu^{(k)}}}      
\newcommand{\xk}{{x^{(k)}}}
\newcommand{\nuk}{{\nu^(k)}}  
\newcommand{\ckE}{{\check E}}
\newcommand{\ckW}{{\check W}}
\newcommand{\hE}{{\hat E}}   
\newcommand{\hH}{{\hat H}}
\newcommand{\hD}{{\hat D}} 
\newcommand{\hDh}{{{\hat D}^{1/2}}}
\newcommand{\hDhn}{{{\hat D}^{-1/2}}}
\newcommand{\sgn}{{\rm sgn}}
\newcommand{\tv}{{\tilde v}}
\newcommand{\te}{{\tilde e}}
%
\journalname{Mathematical Programming}
%
\begin{document}

\title{A modified nearly exact method for solving low-rank trust region subproblem%
\thanks{This work was partially supported by
NSF Grants CCR-0203113, INT-9910084 and CCF-0430644, and 
ONR Grants N00014-03-1-0401 and N00014-05-1-0183.}
}
\titlerunning{A modified nearly exact method for solving LRTR subproblem} % if too long for running head
\author{Zhaosong Lu \and 
        Renato D. C. Monteiro%
}

\institute{Zhaosong Lu \at
              Department of Mathematical Sciences, Carnegie Mellon University, 
              Pittsburgh, Pennsylvania 15213-3890, USA \\
              \email{zhaolu@andrew.cmu.edu}           %  \\
           \and
           Renato D. C. Monteiro  \at
              School of ISyE, Georgia Institute of Technology, 
              Atlanta, Georgia 30332, USA 
              \email{monteiro@isye.gatech.edu}
}

\date{June 1, 2004 \\
1st revision: December 8, 2004 \\
2nd revision: September 16, 2005}
% The correct dates will be entered by the editor

\maketitle

\begin{abstract}
In this paper, we first discuss how the nearly exact (NE) method proposed by
Mor{\'e} and Sorensen \cite{MORE83} for solving trust region subproblems
can be modified to solve large-scale ``low-rank'' trust region TR subproblems
efficiently. Our modified algorithm completely avoids computation of
Cholesky factorizations by instead relying primarily on the Sherman-Morrison-Woodbury
formula for computing inverses of ``diagonal plus low-rank'' type matrices.
We also implement a specific version of the modified log-barrier (MLB)
algorithm proposed by Polyak \cite{POLYAK92} 
where the generated log-barrier subproblems are solved by a trust region method. 
The corresponding direction finding TR subproblems are of the low-rank
type and are then solved by our modified NE method.
We finally discuss the computational results of our implementation
of the MLB method and its comparison with a version of LANCELOT \cite{CONN92}
based on a collection extracted from CUTEr \cite{GoOrTo03}
of nonlinear programming problems with simple bound constraints.
\keywords{nearly exact method \and trust region method \and large-scale optimization 
\and limited-memory BFGS method \and Sherman-Morrison-Woodbury formula}
\subclass{90C06 \and 90C30 \and 65K05}
\end{abstract}

\section{Introduction}
\label{introduction}
Trust region algorithms are classical methods
for solving both convex and nonconvex nonlinear optimization problems.
They are known to possess strong convergence properties
(see Fletcher \cite{FLETCHER80}). At each iteration of a trust region
method, the following is specified:
i) a ``simple'' approximation $\phi(\tilde x)$ to the objective function,
called the model, and; ii) a region $T$ around
the current iterate $x$, where $\phi(\tilde x)$ is believed to provide
a good approximation to the objective function.
An approximate solution $p$ of the subproblem
$\min_p \{ \phi(x+p) : x+p \in T\}$ is then computed,
and the next iterate $\hat x$ is set to be $\hat x:=x+p$ provided
there is a ``significant'' objective function progress;
otherwise, we define the next iterate $\hat x$ as $\hat x :=x$. 
In both cases, the region $T$ might be updated and the process 
is then repeated until a desirable iterate is obtained.

In most trust region methods, the above subproblem is either of
or reduces to the following form:
\beq \label{subproblem}
\begin{array}{l}
\mbox{minimize} \,\{ \,q(p): \|p\|_M \leq \Delta \,\} 
\end{array}
\eeq
where $\Delta$ is a positive parameter, $M$ is a symmetric positive-definite 
matrix referred to as the scaling matrix, $\|\cdot\|_M$ is the $M$-$norm$ 
defined as
\[
\|x\|_M = \sqrt{x^T M x}, \ \ \ \forall x \in \Re^n, 
\]
and $q: \Re^n \to \Re$ is the quadratic function defined as  
\beq \label{model}
q(p) = g^T p + \frac{1}{2} p^T H p, \ \ \ \forall p \in \Re^n,    
\eeq
for some $g \in \Re^n$ and symmetric matrix $H \in \Re^{n \times n}$.
The matrix $H$ can be either the 
Hessian of the objective function or some approximation of it.
%Recently a deep study on the TR subproblem \eqnok{subproblem} is provided
%in Fortin and Wolkowicz \cite{FORTIN02}. It is extended to semidefinite 
%programming(SDP) framework by using modern 
%primal-dual approaches. It is proved to be always a well-posed 
%problem. 
%Some efficient methods are given to solve it under SDP framework.        

There are at least three well-known methods available in the literature for finding 
an ``approximate'' solution of TR subproblem \eqnok{subproblem}, which achieves at least as 
much reduction in the model $q$ as the reduction achieved by the so called ``Cauchy point'' 
(see for example Mor{\' e} \cite{MORE831} and Chapter 4 of Nocedal and Wright 
\cite{NOCEDAL99}). The first method is the {\it dogleg} method proposed by Powell
 \cite{POWELL70}, and later modified by Dennis and Mei \cite{DeMe79}, which is 
appropriate when the model Hessian $H$ is positive definite. Recently, Zhang and 
Xu \cite{ZhXu99} proposed a {\it dogleg} method for the case when $H$ is indefinite, 
which is based on the estimation of the most negative eigenvalue of $H$ and computation 
of the stable Bunch-Parlett factorization of $H$ (see \cite{BunPar71}). The second method 
is the {\it two-dimensional subspace minimization} method proposed by Shultz et al.\ 
\cite{ShScBy85}, which can be applied when $H$ is indefinite, though it also 
requires an estimate of the most negative eigenvalue of $H$. 
%The third method, due to Steihaug \cite{STEIHAUG83}, is most appropriate when $H$ is 
%the Hessian of the objective function and when this matrix is large and sparse.
The third method is the {\it Steihaug} method proposed independently 
by Steihaug \cite{STEIHAUG83} and Toint \cite{Toint81}, which is most appropriate 
when $H$ is the Hessian of the objective function and when this matrix is large and 
sparse. Based on similar ideas as in the {\it Steihaug} method,
Gould et al.\ \cite{GoLuRoTo99} proposed a method
whose resulting TR subproblems are tridiagonal and
can be efficiently solved by the
Mor{\' e} and Sorensen algorithm \cite{MORE83}.

Besides the three ``approximate'' methods mentioned above, 
there is a method due to Mor{\'e} and Sorensen \cite{MORE83}, 
which finds an approximate solution of the TR
subproblem \eqnok{subproblem} in a stronger 
sense (see \eqnok{approximate}).
Following standard convention, we will refer to such  
solutions as ``nearly exact'' (NE) solutions
and to methods for computing them as 
NE methods. Since the NE method of \cite{MORE83} 
requires repeated computations of Cholesky factorizations of diagonal 
displacements of $H$, it is suitable only for small- to medium-sized problems. 
%However, if the matrices $H$ and $M$ have the 
%special structure as described in \eqnok{hessian} and \eqnok{scale} below, 
%we will show that the nearly exact method can be well modified such that 
%the resulting method is efficient and robust for the large-scale problems. 
%This is the main goal of this paper.
 
The main goal of this paper is to develop a method for computing
NE solutions of the TR subproblem \eqnok{subproblem}
when $H$ and $M$ are large-scale matrices having the following
special structures:
%In this paper, we will consider the specialized TR subproblem, where the matrices 
%$H$ and $M$ are the positive diagonal matrix plus explicit low-rank
%update, namely,
\beqa
H &=& D + V E V^T, \label{hessian} \\
M &=& \tD + \tV \tE {\tV}^T \succ 0, \label{scale}
\eeqa
where $D$, $\tD$ and $\tE$ are positive diagonal matrices, $V$ and 
$\tV$ have few number of columns (say less than 10), and $E$ 
is a diagonal matrix. We will refer to the resulting subproblem as
the ``low-rank trust region'' (LRTR) subproblem.
We will show that every step of the NE method of \cite{MORE83}
can be properly modified to handle the LRTR subproblem and
also that the resulting modified NE (MNE) method is quite efficient and robust
for computing NE solutions of large-scale LRTR subproblems.

LRTR subproblems arise in several contexts.  
For example, when using trust region methods to solve 
unconstrained or linear-equality constrained
minimization problems, the matrix $H$ is usually
obtained by using a low-rank update (memoryless) formula and
the resulting $H$ has the structure specified
in \eqnok{hessian}. In such a case, the scaling matrix $M$ is chosen
as either the identity matrix or some other positive definite matrix
whose structure is as specified in \eqnok{scale} and
depends on the specific problem at hand.
It is well known that many constrained minimization problems can be solved by 
minimizing a sequence of unconstrained ones, obtained by using either
the penalty, log-barrier, augmented Lagrangian multiplier
(see for example \cite{NOCEDAL99}), or modified log-barrier methods
(see Polyak \cite{POLYAK92}). Thus, the NE method developed in this 
paper for solving the LRTR subproblem can potentially be used in solving 
many optimization problems.

%An alternative TR method for large-scale unconstrained nonlinear 
%programming is the {\it subspace trust region method}, independently 
%proposed by Wang, Wen and Yuan \cite{WaWeYu04} (see also \cite{Yuan04}). 
%Their method is based on a subspace trust region method and limited memory 
%BFGS method. 

Independently to the present work, Wang, Wen and Yuan \cite{WaWeYu04} 
(see also \cite{Yuan04}) proposed the {\it subspace trust region method} 
for large-scale unconstrained nonlinear optimization problems, which is 
based on a subspace trust region method and limited memory BFGS method. 
The size of the TR subproblems in their method is very small 
(usually about 6 to 16). These TR subproblems can be solved very 
quickly and exactly, and hence their methods can be used to 
solve large-scale unconstrained nonlinear problems.
The method proposed in \cite{WaWeYu04} is generally different from our limited 
memory quasi-Newton TR method based on solving a sequence of 
full-dimensional LRTR subproblems since the TR subproblem of one 
method generally cannot be converted into an equivalent one of
the other method.
%Also, when the hard case (see pp.\ 6) occurs in our method, the solution
%of our LRTR subproblem does not necessarily lie in the low dimensional
%subspace described in equation (1.2) of Wang, Wen and Yuan 
%\cite{WaWeYu04}.

The following notations are used throughout our paper.
We denote the $k$-th coordinate vector by $e_k$ and the identity matrix 
by $I$. Their dimensions are always clear from the context.
The symbol $\Re^n$ denote the $n$-dimensional Euclidean space.
The set of all $m \times n$ matrices with real entries
is denoted by $\Re^{m \times n}$.
For $J \subseteq \{1,\ldots,n\}$ and $w \in \Re^n$, we let $w_J$ 
denotes the subvector $[w_i]_{i \in J}$;
moreover, if $E$ is an $m \times n$ matrix then $E_J$ denotes
the $m \times |J|$ submatrix of $E$ corresponding to $J$.
For a vector $w \in \Re^n$, $\Diag(w)$ denote the diagonal matrix 
whose $i$-th diagonal element is $w_i$ for $i=1,\ldots,n$, and for 
any real number $\alpha$, $w^{\alpha}$ denote the vector whose 
$i$-th component is $w^{\alpha}_i$ whenever it is well-defined 
for $i=1,\ldots,n$. The Euclidean norm, the $1$-norm and the 
$\infty$-norm are denoted by $\|\cdot\|$, $\|\cdot\|_1$ and
$\|\cdot\|_{\infty}$, respectively. For a matrix $E$, 
$\R(E)$ denotes the subspace generated by the columns of $E$ and 
$\Ker(E)$ denotes the subspace orthogonal to the rows of $E$. The
superscript $^T$ denotes transpose. For any real
symmetric matrix $E$, $\lambda_{\min}(E)$ (resp., $\lambda_{\max}(E)$) 
denotes the minimal (resp., maximal) eigenvalue of the matrix $E$; 
$E \succeq 0$ (resp., $E \succ 0$) denotes that $E$ is positive 
semi-definite (resp., positive definite).

Before ending this section, we provide one example to show how
the LRTR subproblem naturally arises in the context of solving
linearly constrained minimization problems using a log-barrier
approach. Indeed, consider the problem
\beq \label{example1}
\begin{array}{ll}
\mbox{minimize} & f(x) \\ [6pt]
\mbox{subject to} & Ax=b, \\
& l \le x \le u,
\end{array}
\eeq
where $f(x)$ is twice continuously differentiable in $\Re^n$, 
$A \in \Re^{m \times n}$ has full row rank, $l, \, u \in \Re^n$ may 
have some components equal to $-\infty$ or $+\infty$, and $m$ is
small.
%Assume that an initial strictly feasible point $x^0$ is available. 
The log-barrier approach applied to \eqnok{example1} consists
of solving the following sequence of log-barrier subproblems
parametrized by $\mu >0$:
\beq \label{example1a}
\begin{array}{ll}
\mbox{minimize} & \phi_\mu(x) \\ [6pt]
\mbox{subject to} & Ax=b,
\end{array}
\eeq              
where 
\[
\phi_\mu(x) = f(x) - \mu \sum_{i=1}^{n} \log(x_i - l_i) -  \mu 
\sum_{i=1}^{n} \log(u_i - x_i).  
\]
Assume that $x$ denotes the current iterate towards a (local)
solution of \eqnok{example1a}. To find the next iterate,
a typical TR method in this context
computes the (potential) displacement $p$ by solving the TR subproblem
\beq \label{example1b}
\begin{array}{ll}
\mbox{minimize} & g^T p + \frac{1}{2} p^T H p \\ [6pt]
\mbox{subject to} & Ap=0,  \\ [3pt]
& \|W^{-1} p\| \le \Delta,
\end{array}
\eeq  
where $W = \Diag([(x-l)^{-2} +(u - x)^{-2}]^{-1/2}) \succ 0$, 
$g =\nabla \phi_\mu(x)$, $\Delta > 0$, and
$H=D+ V E V^T$ is an approximation to $\nabla^2 \phi_\mu(x)$ obtained
by using a low-rank (memoryless) update formula.
Thus, $H$ has the structure as in \eqnok{hessian}, 
%where the matrix $D$ is positive diagonal, 
%the matrix $V$ has a few number of columns, and the matrix $E$ is diagonal. 
%Indeed, $H$ has the same structure as described in \eqnok{hessian}, 
Now, let $B \subset \{1,\ldots,n\}$ be a basic index set and let $N$
denote its complement. By permuting the columns of $A$, $W$ and $D$,
we may assume that
%$N=\{1,\ldots,n-m\}$ and 
%$B = \{n-m+1, \ldots, n\}$. Then we have
\[
A = [A_B, \, A_N], \ \  W = \Diag(W_B,\, W_N), \ \ D = \Diag(D_B, D_N).
\]
%where $W_B, \, W_N, \, D_B,\, D_N$ are the sub-diagonal matrices of 
%$W$ and $D$ corresponding to the index sets $B$ and $N$, respectively. 
Since $Ap = 0$, we have $p_B = -A_B^{-1}A_N p_N$. Eliminating $p_B$ from
\eqnok{example1b}, we obtain the following equivalent LRTR subproblem
\beq \label{example1c}
%\mbox{(LRTR)} \qquad 
\begin{array}{ll}
\mbox{minimize} &  {\bar g}^T p_N + \frac{1}{2} p^T_N 
\bar H p_N \\ [6pt]
\mbox{subject to} & \|p_N\|_M \le \Delta,
\end{array}
\eeq       
where 
\[
\bar g = S^T g, \ \ \ \bar H = S^T H S, \ \ \ M = S^T W^{-2} S, \ \ \ S =
\left[\begin{array}{c}
      -A_B^{-1}A_N \\ 
       I  \end{array} \right].
\]
Thus, we easily see that
\beqas
M 
%&=& S^T W^{-2} S, \nn \\
%&=& [I -(A_B^{-1}A_N)^T] \Diag(W_N^{-2}, W_B^{-2}){\left[\begin{array}{c}
%     I \nn \\
%     -A_B^{-1}A_N \end{array} \right]}, \nn \\ 
&=& W_N^{-2} + (A_B^{-1}A_N)^T W_B^{-2} (A_B^{-1}A_N), \\ [6pt]
\bar H 
%  &=& S^T H S, \nn \\ [6pt]
%  &=& S^T (D+VEV^T)S, \nn \\
%  &=& [I -(A_B^{-1}A_N)^T] \Diag(D_N, D_B) {\left[\begin{array}{c}
%       I \nn \\
%       -A_B^{-1}A_N \end{array} \right]} 
%       +(S^TV)E(V^TS), \nn \\
&=& D_N +(A_B^{-1}A_N)^T D_B (A_B^{-1}A_N) +(S^TV)E(V^TS).  
\eeqas
Noting that $A$ has low full rank and $H$ has the structure specified
in \eqnok{hessian}, we immediately see that the matrices
$M$ and $\bar H$ themselves have the structure as 
in \eqnok{scale} and \eqnok{hessian}, respectively.
Thus, the subproblem \eqnok{example1c} is indeed an LRTR subproblem.
%Actually, a lot of unconstrained 
%problems are the special cases of \eqnok{example1}. 
%Moreover, many constrained problems can be reformulated as 
%unconstrained ones. Therefore, the modified nearly exact method 
%developed in this paper are applicable for many optimization 
%problems.
         
The outline of the paper is as follows. In Section \ref{2}, we 
review the NE method proposed by Mor{\'e} and Sorensen 
\cite{MORE83}. In Section \ref{3}, we discuss how this method can be 
modified in order to solve large-scale LRTR subproblems efficiently.
In Section \ref{4}, we first review the modified 
log-barrier (MLB) algorithm proposed by Polyak \cite{POLYAK92} and 
implement a specific version of this algorithm where the generated
log-barrier subproblems are solved by a trust region method whose
direction finding subproblems are of the LRTR type.
The LRTR subproblems are then solved by our modified NE method.
Section \ref{4} also gives computational results of our implementation
of the MLB method and its comparison with a version of LANCELOT \cite{CONN92}
based on a collection extracted from CUTEr \cite{GoOrTo03}
of nonlinear programming problems with simple bound constraints.

\section{Review the NE method for solving TR subproblem}
\label{2}

%As we know, the performance of trust region method strongly depends on how to 
%solve TR subproblem \eqnok{subproblem}. On one hand, finding the exact solution of
%\eqnok{subproblem} may be very time consuming. On the other hand,
It is well-known that the TR subproblems which arise in a TR method
does not need to be solved exactly to guarantee the global
convergence of the algorithm. For example, it has been shown by Mor{\'e}   
and Sorensen \cite{MORE83} (see also \cite{MORE831}) that,
under some mild conditions, good theoretical and numerical
convergence results for a standard TR method can be obtained if
$p$ is chosen so that
\beq \label{approximate}
q(p) \le \tau_1 q^* \ \ \ \mbox{and} \ \ \ \|p\|_M \le \tau_2 \Delta
\eeq
for some positive constants $\tau_1$ and $\tau_2$, 
where $q^*$ is the optimal 
value of TR subproblem \eqnok{subproblem}.
(Note that $q^* \le 0$ and that $q^*=0$ if and only if
$g = 0$ and $H \succeq 0$.)
We will refer to such vectors $p$ as NE solutions of
\eqnok{subproblem}.

The NE method proposed by Mor{\'e} and Sorensen \cite{MORE83} 
is a method for computing a NE solution $p$ of \eqnok{subproblem}. 
In this section, we review the technical 
results of the NE method of \cite{MORE83} for solving TR
subproblem \eqnok{subproblem}
(see \cite{CONN00,FORTIN02,MORE83,MORE831,NOCEDAL99} for more details).
The computational difficulties of using the NE method to
TR subproblems corresponding to large-scale optimization problems
are also presented. But, in Section \ref{3}, we will show that
the NE method can be suitably modified to overcome these difficulties
if all the TR subproblems are constructed as LRTR ones.

This section is divided into five subsections. In Subsection 
\ref{2.1}, we discuss the necessary and sufficient optimality 
conditions for a global solution of the TR subproblem 
\eqnok{subproblem}. In Subsection \ref{2.2}, we discuss some 
classical and easily verifiable sufficient conditions for $p$ 
to be a NE solution of \eqnok{subproblem} (see for example 
\cite{CONN00} and \cite{MORE83}). In Subsection \ref{2.3}, 
we discuss how Newton method applied to a classical one dimensional 
nonlinear equation provides an estimate of the optimal Lagrange 
multiplier associated with the constraint of \eqnok{subproblem}.
Since the search for the optimal Lagrange multiplier requires
the estimation of ever-improving lower bounds for it, we discuss
in Subsection \ref{2.4} how these bounds are normally generated.
The complete NE method of \cite{MORE83} and its computational 
difficulties in the context of large-scale problems are also 
discussed in this subsection.

\subsection{Characterization of the solution of TR subproblem}
\label{2.1}

In this subsection, we provide optimality conditions
which characterize the global
solutions of subproblem \eqnok{subproblem}.

The proof of the next lemma, which provides the above mentioned
optimality conditions, is given in Theorem 7.4.1 on pp.~201 of 
Conn et al.\ \cite{CONN00}. 
(This result was obtained independently by 
Gay \cite{GAY81} and Sorensen \cite{SORENSEN82a}.)

\begin{lemma} \label{optcond}
$p$ is a global solution of TR subproblem \eqnok{subproblem} if only 
if $\|p\|_M \le \Delta$ and there exists $\lambda \ge 0$ such that
\beqa 
H(\lambda) p &=& -g, \label{equation} \\
\lambda(\Delta - \|p\|_M) &=& 0, \label{complement} \\
H(\lambda) &\succeq& 0, \label{posdef}     
\eeqa   
where $H(\lambda) \equiv H + \lambda M$ for any $\lambda \in \Re$. 
Moreover, there exists a unique $\lambda^* \ge 0$ such that:
\begin{itemize}
\item[i)]
$\lambda=\lambda^*$ for every pair $(p,\lambda)$ as above;
\item[ii)]
if $H(\lambda^*) \succ 0$ then \eqnok{subproblem} 
has a unique global optimal solution.  
\end{itemize}
\end{lemma}

We now introduce some notation. Define
\beq \label{lambda1}
\lambda_1 \equiv \lambda_{\min}(\Mhn H \Mhn),  \ \ \ \
\hat {\lambda} \equiv \max(-\lambda_1, 0).
\eeq
%Moreover, we define the critical vector $p_{\rm crt}$ as
%\[
%p_{\rm crt} \equiv -H(\hat \lambda)^\dag g,
%\]
%where the superscript $^\dag$ denotes the Moore-Penrose
%generalized inverse. It can be shown that
%$p_{\rm crt} = \lim_{\lambda \downarrow  \hat \lambda} p(\lambda)$,
Moreover, we define
\beq \label{step1}
p(\lambda) \equiv -H(\lambda)^{-1} g,
\eeq
for every $\lambda \in \Re$ for which the above inverse exists. 
It is well-known that, when $g \ne 0$, the function
$\|p(\lambda)\|_M$ is strictly decreasing and convex on 
$(\hat \lambda, \infty)$.    

We now describe how an exact solution for
\eqnok{subproblem} can be computed, depending on which of
the following three cases occur:
\begin{itemize}
%\item[{\bf Case 1}]
\item[1)]
If there exists $\tilde \lambda \in (\hat \lambda, +\infty)$ such 
that $\tilde \lambda$ solves the equation
\beq \label{newton1}
\|p(\lambda)\|_M = \Delta,
\eeq
then $\lambda^*=\tilde \lambda >0$
and $p(\tilde \lambda)$ is the unique solution 
of \eqnok{subproblem}. This case is usually referred to as 
the ``easy'' one. (Note that this case occurs if and only 
if $\lim_{\lambda \downarrow \hat \lambda}\|p(\lambda)\|_M > \Delta$.)
\item[2)]
If $\lambda_1 >0$ and \eqnok{newton1} has no solution in 
$(\hat \lambda, +\infty)$, then $\lambda^* = \hat \lambda =0$ and
$p(0)$ is a solution of \eqnok{subproblem}. (This case
can easily be detected and usually referred to as
the ``interior convergence'' one.) 
\item[3)]
If $\lambda_1 \le 0$ and \eqnok{newton1} has no solution in
$(\hat \lambda, +\infty)$, then $\lambda^*=\hat \lambda=-\lambda_1$.
Hence, there exist $0 \neq u \in \Ker(H+\hat \lambda M)$ and 
$\alpha^M \in \Re$ such that 
\[
\|p_{\rm crt} + \alpha^M u\|_M = \Delta,
\]
where $p_{\rm crt} \equiv -H(\hat \lambda)^\dag g$, and the superscript 
$^\dag$ denotes the Moore-Penrose generalized inverse.
Using Lemma \ref{optcond}, one easily sees that
$p_{\rm crt} + \alpha^M u$ is a solution of \eqnok{subproblem}.
This case is usually referred to as the ``hard'' one.
\end{itemize}     

Some steps in the NE algorithm for solving \eqnok{subproblem}
require to test whether $\lambda > \hat {\lambda}$ or
$\lambda > \lambda^*$. These two inequalities can be checked by using
the following easily verifiable characterizations:
i) $\lambda > \hat {\lambda}$ if and only if
$\lambda >0 $ and $H(\lambda) \succ 0$; ii)
$\lambda > \lambda^*$
if and only if $\lambda > \hat {\lambda}$ and
$\|p(\lambda)\|_M < \Delta$.

\subsection{Termination conditions}
\label{2.2}

Among the three cases mentioned in Subsection \ref{2.1},
only the interior convergence case can be implemented exactly.
When the other two cases occur, we can only expect to
obtain an approximate solution of \eqnok{subproblem}.
In this subsection, we review some sufficient conditions
for a vector $p \in \Re^n$ to be a NE solution of \eqnok{subproblem}
when either the easy or hard case occurs.

While looking for a scalar $\lambda>\hat \lambda$ satisfying
\eqnok{newton1}, we might simply stop when
$\left|\,\|p(\lambda)\|_M -\Delta \,\right| \le k_e \Delta$, where 
$k_e \in (0,1)$ is a fixed tolerance. In this case, the following result
establishes that $p(\lambda)$ is a NE solution of \eqnok{subproblem}.
Its proof is similar to the one given in Lemma 7.3.5 on
pp.~195 of \cite{CONN00} (see also \cite{MORE83}).

\begin{lemma} \label{easy}
If $\lambda > \hat \lambda$ satisfies
$\left|\,\|p(\lambda)\|_M -\Delta \, \right| \leq k_e \Delta$
for some $k_e \in (0,\, 1)$, then
$q(p(\lambda)) \leq (1- k_e)^2 q^*$.
\end{lemma}

The following result describes how an approximate version of the
hard case yields NE solutions for \eqnok{subproblem}.
Its proof is similar to the one given in Lemma
7.3.6 on pp.~196 of \cite{CONN00} (see also \cite{MORE83}).

\begin{lemma} \label{hard}
Suppose that $\lambda > \lambda^*$, $\alpha \in \Re$ and 
$u \in \Re^n$ such that $\|u\|_M =1$ satisfy
\beq \label{hardterm}
\alpha^2 u^T H(\lambda) u \leq k_h \left(p(\lambda)^T
H(\lambda) p(\lambda) + \lambda \Delta^2 \right) \,\, \mbox{and} \,\,
\|p(\lambda) + \alpha u\|_M = \Delta,
\eeq
for some $k_h \in (0,\, 1)$. Then,
$q(p(\lambda)+ \alpha u) \leq (1- k_h) q^*$.
\end{lemma}

The main use of Lemma \ref{hard} is related with the hard case, but we 
emphasize that the result can also be applied in other cases. We 
now explain how to satisfy the conditions of Lemma \ref{hard} in the hard 
case. From the discussion in Subsection \ref{2.1}, we know that 
$\lambda^*=-\lambda_1$. Hence, if $\lambda-\lambda^*$ is sufficiently 
small, then $H(\lambda)$ is nearly singular and thus there exists 
a vector $u$ such that $\|u\|_M = 1$ and $u^T H(\lambda) u$ 
is nearly zero. Once the pair $(\lambda,u)$ is determined,
a scalar $\alpha$ satisfying the equality in \eqnok{hardterm}
can be easily obtained by solving the following problem
\beq \label{qp1}
\begin{array}{ll}
\mbox{minimize}_{\,\alpha} & \frac{1}{2} (p+\alpha u)^T H (p+\alpha u)
+ g^T(p+\alpha u)\\ [6pt]
\mbox{subject to} & \|p+\alpha u\|_M = \Delta,
\end{array}
\eeq
where $p=p(\lambda)$.
%We next provide a formula for computing such a $\alpha$. Let $\bar p = \Mh p$, 
%\ $\bar u = \Mh u$, \ $\bar H = \Mhn H \Mhn$ and   
%$\bar g = \Mhn g$. Then the problem \eqnok{qp1} is equivalent to
%\beq \label{qp2}
%\begin{array}{ll}
%\mbox{minimize}_{\,\alpha} & \frac{1}{2} (\bar p+\alpha \bar u)^T
%\bar H (\bar p+\alpha \bar u) + {\bar g}^T(\bar p+\alpha \bar u)\\ [6pt]
%\mbox{subject to} & \|\bar p+\alpha \bar u\| = \Delta.
%\end{array}
%\eeq
%Using the result of the optimal solution of \eqnok{qp2} from \cite{MORE83}, 
%\[
%\alpha = \frac{\Delta^2 -\|\bar p\|^2}{{\bar p}^T \bar u
%+ \sgn ({\bar p}^T \bar u)[({\bar p}^T \bar u)^2+ \Delta^2
%-\|\bar p\|^2]^{1/2}},
%\]
%together with $\bar p = \Mh p$ and $\bar u = \Mh u$, 
Using the well-known formula of $\alpha$ for the case when $M=I$
(see page 558 of \cite{MORE83}), it is easy to see that
the optimal solution of \eqnok{qp1} is given by
\beq \label{alpha}
\alpha = \frac{\Delta^2 -\|p\|^2_M}{p^TM u + \sgn(p^T M u)[(p^T M u)^2
+ \Delta^2 -\|p\|^2_M]^{1/2}},
\eeq
where the function $\sgn: \Re \to \Re$ is defined as $\sgn(t)=1$ if 
$t \geq 0$, and $\sgn(t)=-1$ if $t<0$.
It can be easily verified that the right hand side
of the inequality in \eqnok{hardterm} evaluated at
a triple $(\lambda,u,\alpha)$ obtained as above stays bounded
away from zero as $\lambda-\lambda^*$ approaches zero.
%Moreover, using \eqnok{step1}, we have
%\[
%p(\lambda)^T H(\lambda) p(\lambda)=g^T H^{-1}(\lambda)g \ge 
%\frac{\|g\|^2_{M^{-1}}}{\lambda+\lambda_n},
%\]
%where $\lambda_n=\lambda_{\max}(\Mhn H \Mhn)$. Hence, 
%if $\lambda-\lambda^*$ is sufficiently small, then  
%$p(\lambda)^T H(\lambda) p(\lambda) > \delta$ for some constant 
%$\delta >0$. From $\|p(\lambda) + \alpha u\|_M = \Delta$, we also 
%see that $\alpha$ is bounded if $\lambda-\lambda^*$ is sufficiently 
%small.
Thus, as $\lambda-\lambda^*$ approaches zero, a triple
$(\lambda,u,\alpha)$ obtained as above will eventually satisfy
\eqnok{hardterm} when the hard case occurs.

The key part to find a triple $(\lambda,u,\alpha)$ satisfying
\eqnok{hardterm} in the hard case is the computation of a
vector $u=u_\lambda$ such that
$u^TH(\lambda)u$ approaches zero as $\lambda \downarrow \lambda^*$.
The NE method of \cite{MORE83} computes such a vector $u$
by first computing the Cholesky factor of $H(\lambda)$, and then
using the LINPACK technique \cite{ClMoStWi79}
(see also Appendix of \cite{MORE83}) for estimating 
its smallest singular value. On the other hand, since our approach
for solving the low rank version of \eqnok{subproblem} does not rely
on computation of Cholesky factorizations, it uses an entirely different
approach to compute a vector $u$ as above (see
Subsection \ref{3.3}).

\subsection{Newton update for $\lambda$}
\label{2.3}

We have seen in the Subsection \ref{2.1} that
$\lambda^*$ is a root of \eqnok{newton1} in the easy case.
Hence, given an approximation $\lambda > -\lambda_1$
of $\lambda^*$, it is natural to try to perform a Newton
iteration at $\lambda$ with respect to \eqnok{newton1}
to obtain a new approximation of $\lambda^*$.
In this subsection, we describe the details of
a Newton iteration applied to a reformulation
of the nonlinear equation \eqnok{newton1} and
discuss its main properties.
 
Since the function $\|p(\lambda)\|_M$ goes to infinity as
$\lambda$ tends to $-\lambda_1$, it is highly nonlinear
near $-\lambda_1$ (see \cite{MORE831} and \cite{NOCEDAL99}).
As a result, Newton method applied directly to
\eqnok{newton1} might not work well when $\lambda$ is near
$-\lambda_1$.  Reinsch \cite{Reinsch71} and
Hebden \cite{Hebden73} independently observed that
Newton method applied to the following alternative reformulation
of \eqnok{newton1} works better in practice:
\beq \label{newton}
\phi(\lambda) \equiv \frac{1}{\|p(\lambda)\|_M} - \frac{1}{\Delta} = 0.
\eeq

The following result describes some important properties of the function 
$\phi(\lambda)$ and provides the formula of a Newton iteration for
\eqnok{newton}.

\begin{lemma} \label{funprop}
Suppose $g \ne 0$. Then, $\phi(\lambda)$ is strictly 
increasing and concave on $(-\lambda_1,\infty)$. Moreover, the Newton 
iterate at $\lambda$ with respect to \eqnok{newton} is 
\beq \label{newtonstep1}
\lambda^+ = \lambda + \left(\frac{\|p(\lambda)\|_M - \Delta}{\Delta}\right)
\left(\frac{\|p(\lambda)\|^2_M} {p(\lambda)^T MH(\lambda)^{-1}M p(\lambda)}\right).
\eeq
\end{lemma} 

\begin{proof}
For $\lambda > -\lambda_1$, let
$\bar p (\lambda) \equiv - \bar H (\lambda)^{-1} \bar g$,
where $\bar g \equiv \Mhn g$ and $\bar H (\lambda) \equiv \Mhn H(\lambda) \Mhn =
\Mhn H \Mhn + \lambda I$. By \eqnok{step1} we have
$\bar p(\lambda) = \Mh p(\lambda)$, which together with \eqnok{newton}
implies that
\beqas
\phi(\lambda) = \frac{1}{\|\bar p(\lambda)\|} - \frac{1}{\Delta}.
\eeqas
Hence, by Lemma 7.3.1 on pp.~183 of \cite{CONN00}, it follows that
the function $\phi(\lambda)$ is strictly increasing and
concave when $\lambda \in (-\lambda_1,\infty)$, and that its
first derivative is given by
\beq \label{deriv2}
\phi'(\lambda) = \frac{{\bar p(\lambda)}^T {\bar H(\lambda)}^{-1} 
\bar p(\lambda)}{\|\bar p(\lambda)\|^3}.
\eeq    
The formula \eqnok{newtonstep1} for the Newton iteration
$\lambda^+=\lambda-\phi'(\lambda)/\phi(\lambda)$ can be easily
derived using \eqnok{deriv2}, $\bar p(\lambda)  = \Mh p(\lambda)$ 
and $\bar H(\lambda) = \Mhn H(\lambda) \Mhn$. 
\qed
\end{proof}  

\vgap

The next result gives a few useful properties of Newton method applied to
\eqnok{newton}.

\begin{proposition} \label{newtonprop}
Suppose $g \ne 0$.
Then the following statements hold:
\bi
\item[a)]
Suppose $\lambda \in (-\lambda_1,\lambda^*)$. Then all Newton
iterates starting from $\lambda$ will stay in
$(-\lambda_1,\lambda^*)$ and converge to the solution
$\lambda^*$ of the equation \eqnok{newton} monotonically. The
convergence is globally Q-linear with the ratio at least
\[
\gamma_{\lambda} = 1 - \frac{\phi'(\lambda^*)}{\phi'(\lambda)} < 1
\]
and is ultimately Q-quadratic.
\item[b)]
Suppose $\lambda \in (\lambda^*,\infty)$. Then the next Newton
iterate $\lambda^+ \in (-\lambda_1,\lambda^*]$ or $\lambda^+
\in (-\infty, -\lambda_1]$.
\ei
\end{proposition}

\begin{proof}
The proof is similar to the ones given in 
Lemmas 7.3.2 and 7.3.3 on pp.~185-186 of \cite{CONN00}. 
\qed
\end{proof}

\vgap

To compute the Newton iterate $\lambda^+$ according to
\eqnok{newtonstep1}, the NE method of \cite{MORE83} 
first computes the lower Cholesky factor $L$ of $H(\lambda)$, and uses
it to first compute a vector $p$ such that
$L L^T p = -g$ and then a vector $w$ such that $Lw = M p$.
By \eqnok{newtonstep1}, we then have
\beq \label{newtonstep2}
\lambda^+ = \lambda + \left(\frac{\|p\|_M-\Delta}{\Delta}\right)
\left( \frac {\|p\|^2_M}{\|w\|^2}\right).
\eeq
In our approach for solving the low rank version of \eqnok{subproblem}
we entirely avoid the computation of Cholesky factorizations by instead
computing the inverse of $H(\lambda)$ by means of the Sherman-Morrison-Woodbury 
(SMW) formula (see Subsection 3.2).

\subsection{A safeguard Newton method}
\label{2.4}

Since a Newton iteration might result in infeasible iterates
$\lambda^+ \le -\lambda_1$, the NE method
of \cite{MORE83} uses some safeguard strategies
to handle such iterates in order to obtain a globally
convergent method for obtaining a NE solution of \eqnok{subproblem}.
The basic idea used is to
bracket $\lambda^*$ by a lower bound $\lambda^L$ and
an upper bound $\lambda^U$ and reduce the length of the interval
$[\lambda^L,\lambda^U]$ by using a clever bisection strategy.
In this subsection, we discuss the details of this hybrid method.

At every iteration of the method, we have two scalars
$\lambda^L$ and $\lambda^U$ such that
$0 \le \lambda^L  \le \lambda^* \le \lambda^U$
and a current approximation $\lambda \in [\lambda^L,\lambda^U]$ of
$\lambda^*$. Each iteration of the method then consists of
updating the quantities $\lambda^L$, $\lambda^U$ and $\lambda$.
We first describe the basic idea used to update $\lambda^L$.
Suppose that $u \in \Re^n$ is a vector such that
$\|u\|_M = 1$. Using \eqnok{lambda1}, we see that for any $\lambda \in \Re$,
\[
u^T H(\lambda) u  = (\Mh u)^T (\Mhn H \Mhn + \lambda I ) (\Mh u)
\ge (\lambda_1 + \lambda) \|u\|_M^2 = \lambda_1 + \lambda.
\]
Defining
\beq \label{lambdaB}
\lambda^B = \lambda^B (\lambda,u) \equiv \lambda - u^T H(\lambda) u,
\eeq
it follows from the above inequality that
$\lambda^B \le - \lambda_1 \le \lambda^*$, or in words,
$\lambda_B$ is a lower bound for $\lambda^*$.
In view of this discussion, we conclude that a natural update for
$\lambdaL$ is simply to let
$\lambdaL \leftarrow \max(\lambdaL,\, \lambda^B)$.

We are now ready to describe how the NE method of \cite{MORE83}
updates the three quantities $\lambda^L$, $\lambda^U$ and $\lambda$.
Fix some constant $\theta \in (0,1)$ (e.g., $\theta = 0.01$).
It is convenient to consider the following three cases
separately:
\begin{itemize}
\item[i)]
Assume $\lambda \le -\lambda_1$. If $\lambda=-\lambda_1$, we perform 
the update $\lambdaL \leftarrow \max(\lambdaL,\, \lambda)$. Otherwise, 
$H(\lambda)$ is indefinite, and hence there exists a vector $u$ such 
that $\|u\|_M=1$ and $u^T H(\lambda)u <0$. One approach to find such 
a vector $u$ is to perform a partial Cholesky factorization of $H(\lambda)$
to find a scalar $\delta > 0$ and a vector $v$ such that
\beq \label{aux}
(H(\lambda)+\delta e_k e_k^T) v = 0 \,\,\, \mbox{and} \,\,\, e_k^T v =1.
\eeq
(for more details, see \cite{MORE83} or pp.~191-192 of \cite{CONN00}).
Letting $u = v/\|v\|_M$ in \eqnok{lambdaB} and using \eqnok{aux},
we easily see that $\lambda^B = \lambda + \delta/\|v\|^2_M$. Then,
we perform the update $\lambdaL \leftarrow \max(\lambdaL,\, \lambda^B)$. 
Finally, we perform the update $\lambda \leftarrow \max(\sqrt{\lambdaL \lambdaU},
\lambdaL+\theta(\lambdaU-\lambdaL))$.
\item[ii)] 
Assume $\lambda \in (-\lambda_1, \lambda^*)$. In this case,
we perform the updates
$\lambdaL \leftarrow \lambda$ and
$\lambda \leftarrow \lambda^+$,
where $\lambda^+$ denotes the Newton iterate defined in 
\eqnok{newtonstep1}.
\item[iii)]
Assume $\lambda > \lambda^*$. In this case, we have seen
in the paragraph following
Lemma \ref{hard} that a vector $u \in \Re^n$ such that $\|u\|_M=1$
can be computed using the LINPACK technique
which makes $u^TH(\lambda)u$ small as long as
$\lambda+\lambda_1$ is small.
This vector $u$, together with $\lambda$, is then used to compute
$\lambda^B$ according to \eqnok{lambdaB}, and the updates
$\lambdaU \leftarrow \lambda$, 
$\lambdaL \leftarrow \max(\lambdaL,\lambda^B)$ and
$\lambda \leftarrow \max(\lambda^+, \lambdaL)$ are then performed.
\end{itemize} 

The above scheme for updating $\lambdaL$, $\lambdaU$ and $\lambda$
can be shown to generate a sequence of $\lambda$'s which approaches
$\lambda^*$. Indeed, if case ii) occurs then the sequence of
$\lambda$'s generated afterwards approaches $\lambda^*$
monotonically from the left in view of Proposition \ref{newtonprop}(a).
Also, Proposition \ref{newtonprop}(b) implies that if case iii) occurs
then either case i) or ii) must occur at the next iteration.
Hence, if case ii) never occurs then case i) must occur infinitely
often. But every time i) occurs, it is easy to see that
the ratio of the length of the interval $[\lambdaL,\lambdaU]$ at the
end of the next iteration and
its length at the current iteration is bounded above
by $\max(\theta, 1-\theta )$. Hence, if case ii) never occurs,
the length of the generated intervals $[\lambdaL,\lambdaU]$ converges
to zero, and thus the generated sequence of $\lambda$'s
approaches $\lambda^*$.

Note that the implementations of cases i) and iii) of the above scheme
for updating $\lambdaL$, $\lambdaU$ and $\lambda$ are based
on the computation of the (partial) Cholesky factorization of a matrix.
In our approach for finding NE solutions of the low rank version of
\eqnok{subproblem}, these cases must be implemented differently
so as to avoid the computation of Cholesky factorizations, which are
known to be expensive for large scale problems.
These alternative implementations of cases i) and iii) are discussed in
detail in Subsections \ref{3.2} and \ref{3.3}.

We are now ready to state the whole algorithm of \cite{MORE83}
for finding a NE solution of \eqnok{subproblem}.

%\begin{minipage}[h]{5in}
\begin{tabbing}
\={\bf Algorithm 1 (NE method for solving \eqnok{subproblem}):}
\\ [4.0pt]
\> \=Let constants $\theta,\,k_e,\,k_h \in (0, 1)$ be given
(e.g., $\theta=0.01,\,k_e=0.1,\, k_h=0.2$). \\ [4pt]
\> \> 1)
\=Find initial scalars $0 \le \lambdaL < \lambdaU$ such that
$\lambda^* \in [\lambdaL,\lambdaU]$ and set 
$\lambda = \lambdaL$. \\ [4pt]
\> \> 2)
Attempt to do Cholesky factorization $H(\lambda) = L L^T$ to check whether 
$\lambda > -\lambda_1$. \\ [3.5pt] 
\> \> 3)
If $\lambda > -\lambda_1$, solve $LL^T p = -g$ for $p$. Check 
for interior convergence and easy termination. \\ [3.5pt]
\> \> \> If $\lambda > \lambda^*$, compute a pair
$(u,\alpha) \in \Re^n \times \Re$ to check for the hard termination
\eqnok{hardterm}. \\[3.5pt]
\> \> 4) 
If $\lambda \le -\lambda_1$, then update $\lambda^L$ and $\lambda$ according 
to case i) as above, and go to step 2). \\ [3.5pt]
\> \> 5)
If $\lambda \in (-\lambda_1, \lambda^*)$, then update $\lambda^L$ and $\lambda$
according to case ii) as above, and go to step 2). \\ [3.5pt]
\> \> 6)
If $\lambda > \lambda^*$, then update $\lambda^L$, $\lambda^U$ and $\lambda$
according to case iii) as above, and go to step 2). \\ [4pt]
{\bf End}
\end{tabbing}
%\end{minipage}

\mgap

For large-scale problems, several computational difficulties arise in
the NE method of \cite{MORE83} above.  
In step 1),  it generally takes ${\cal O} (n^2)$ amount of
arithmetic operations to find the initial $\lambdaU$
(see Section 7.3.8 of \cite{CONN00}), which is somehow 
expensive for large-scale problems. The Cholesky or partial Cholesky 
factorization of $H(\lambda)$ used in steps 2), 4), 6) and 7) needs 
${\cal O} (n^3)$ amount of arithmetic operations. It will be prohibitive 
for large-scale problems. In the next section, 
we will modify the NE method of \cite{MORE83} above to overcome those difficulties 
for solving the LRTR subproblem. 

\section{A modified NE method for solving LRTR subproblem}
\label{3}

In this section a modified NE (MNE) method for solving large-scale LRTR 
subproblems is presented. We follow the framework of Algorithm 1 as described 
in Section \ref{2}. Our main effort is to overcome the computational difficulties 
mentioned in Section \ref{2}. 

This section is divided into four subsections. A more efficient approach for 
checking whether $H(\lambda)$ is positive definite is given in Subsection \ref{3.1}. 
We also modify the approach of solving the linear equation $H(\lambda)p = -g$ and computing 
a Newton iterate in this subsection. A more efficient approach for dealing with hard case  
termination is developed in Subsection \ref{3.2}. In Subsection \ref{3.3}, 
a more efficient approach for improving$\lambdaL$ when $\lambda < \lambda_1$ 
is given. In Subsection \ref{3.4}, we develop a cheaper approach to initialize $\lambdaU$.
We emphasize that all modified approaches completely avoid computing Cholesky 
or partial Cholesky factorization of large-scale matrices. 

\subsection{Checking positive definiteness of $H(\lambda)$ and solving 
$H(\lambda) p = -g$}
\label{3.1}

In step 2) of Algorithm 1, given any $\lambda \geq 0$, the NE method of 
\cite{MORE83} checks whether $H(\lambda)$ is positive definite by computing 
the Cholesky factorization of it, which is very expensive and even prohibitive 
for the large-scale problems. In this subsection, we provide a more efficient 
method instead, which needs ${\cal O}(n)$ amount of arithmetic operation for 
large-scale LRTR problems. Furthermore, in steps 5) and 6) of Algorithm 1, 
the NE method of \cite{MORE83} uses the Cholesky factor of $H(\lambda)$ to 
solve the linear equation $H(\lambda) p = -g$ and compute Newton iterate, 
respectively. We will modify those approaches as well in this subsection. 

The following theorem provides the main tool for the analysis in this subsection.

\begin{theorem} \label{testpsd}
Let $\hE \in {\cal S}^m$, $\hV \in \Re^{n \times m}$ and
an invertible matrix $\hD \in {\cal S}^n$ be given and
define $\hH \equiv \hD + \hV \hE \hV^T$ and $\hW \equiv \hV^T \hD^{-1} \hV$. Then,
the following statements hold:
\begin{itemize} 
\item[i)]
If $\hE$ is invertible, then
$\hH$ is invertible  if and only if $\hE^{-1} + W$ is invertible.
\item[ii)]
If $\hD \succ 0$, then $\hH \succeq 0$ if and only if 
$\hW + \hW \hE \hW \succeq 0$.
\item[iii)]
If $\hD \succ 0$, then $\hH \succ 0$ and $\hV$ has full column rank if
and only if $\hW + \hW \hE \hW \succ 0$.
\end{itemize}
\end{theorem}

\begin{proof}
It is well-known that if $\hE^{-1} + \hW$ is invertible then
the SWM formula applied to $\hH$ implies that
$\hH^{-1} =  \hD^{-1} - \hD^{-1} \hV ( \hE^{-1} + \hW )^{-1} \hV^T \hD^{-1}$.
Similarly, if $\hH=\hD+\hV \hE \hV^T$ is invertible then the SMW formula applied
to $\hE^{-1} + \hW$ reveals that this matrix is invertible. Hence,
i) follows.

To prove statement ii), assume that $\hD \succ 0$.
We can then write $\hH$ as $\hH = \hDh F \hDh$, where
$F \equiv I + \hDhn \hV \hE \hV^T \hDhn$. Clearly, $\hH \succeq 0$
if and only if $F \succeq 0$. In view of the decomposition
$\Re^n = \Ker ( \hV^T \hDhn ) + \R ( \hDhn \hV )$,
it follows that, for any $p \in \Re^n$, there exist $u,\, v \in \Re^n$ and
$y \in \Re^m$ such that
\beqas
p = u + v, \ \ \  \hV^T \hDhn u = 0, \ \ \  v = \hDhn \hV y.
\eeqas
This relation, the facts that $Fu=u$ and $u^Tv=0$,
and some simple algebraic manipulation imply that
\beqas
p^T F p 
&=& (u+v)^T F (u+v) = u^TFu +2u^TFv+v^TFv , \\
&=& \|u\|^2 + v^T F v = \|u\|^2 + y^T (\hW + \hW \hE \hW) y.     
\eeqas 
By the arbitrariness of $p$, we easily see that $F\succeq 0$ 
if only if $\hW + \hW \hE \hW \succeq 0$. This, together with the 
fact that $\hH \succeq 0$ if only if $F \succeq 0$, implies that 
ii) holds. Under additional assumption that $\hV$ has full column 
rank, the statement iii) can be shown by using a similar argument 
as ii). 
\qed 
\end{proof}  

\mgap

We now describe an efficient approach to check whether $H(\lambda)$
with $\lambda \ge 0$ is positive definite in step 2) of Algorithm 1.
Noting that $H(\lambda)=H+\lambda M$, we 
see from \eqnok{hessian} and \eqnok{scale} that
\beq \label{Hlambda}
H(\lambda) = \hD + \hV \hE \hV^T,
\eeq
where 
\beq \label{hDVE}
\hD \equiv D + \lambda \tD \succ 0, \ \ \ \hV \equiv (V, \tV), \ \ \ 
\hE \equiv \Diag(E, \lambda \tE). 
\eeq
Hence, by Theorem \ref{testpsd} ii), if $\hV$ has
full column rank, then we can check the positive
definiteness of $H(\lambda)$ by checking whether
\beq \label{matX1} 
X \equiv \hW+\hW \hat E \hW
\eeq 
is positive definite. 
%where $\hW \equiv \hV^T \hD^{-1} \hV$.
On the other hand, if $\hV$ does not
have full column rank, then we determine a matrix $R$ with full
column rank and a matrix $T$ such that $\hV=R T$.
It then follows that $H(\lambda)=\hD+\hV \hE \hV^T= 
\hD+R \ckE R^T$, where $\ckE \equiv T\hE T^T$.
Hence, Theorem \ref{testpsd} ii) can now be used to check
the positive definiteness of $H(\lambda)$ by checking  
whether $\ckW+\ckW \ckE \ckW \succ 0$, where $\ckW \equiv R^T \hD^{-1} R$.
For convenience of the presentation, we will assume
throughout the remaining subsections that $\hV$ has full column 
rank. 

Recall that one of the requirements for \eqnok{subproblem} to be a LRTR subproblem 
is that the number of columns of $\hV$ be small. In this case, $X$ is a small-sized 
matrix whose positive definiteness can be checked by performing a relatively cheap Cholesky 
factorization. Since the amount of arithmetic operations to compute $X$ for a large-scale 
LRTR subproblem is ${\cal O}(n)$, the above approach for checking whether $H(\lambda)$ is 
positive definite only requires ${\cal O}(n)$ arithmetic operations. If $H(\lambda)$ turns 
out to be positive definite, we can then solve the linear system
$H(\lambda)p = -g$ by means of SMW formula as 
\beqas 
p &=& -{H(\lambda)}^{-1} g = -\left(\hDn - \hDn \hV 
(\hEn + \hV^T \hDn \hV)^{-1} \hV^T \hDn\right)g, \nn \\
&=& - \hDn g + \hDn \hV (\hEn + \hW)^{-1} \hV^T \hDn g. 
\eeqas   
Note that the matrix $\hEn + \hW$ is invertible due to the fact $H(\lambda) \succ 0$ 
and Theorem \ref{testpsd} i). Since the size of $\hEn + \hW$ is small, 
this approach for solving the linear system $H(\lambda)p = -g$ 
also requires ${\cal O}(n)$ arithmetic operations. 
Note that, in the context of a LRTR subproblem, the Newton iterate $\lambda^+$ can be efficiently 
computed by means of \eqnok{newtonstep1}, which requires solving another linear system with 
coefficient matrix $H(\lambda)$. 

%In steps 5) and 6) of Algorithm 1, the NE method of \cite{MORE83} also uses the 
%Cholesky factor of $H(\lambda)$ to compute the Newton iterate (see 
%Subsection \ref{2.3}). We now modify step 4) as follows. Indeed, we first 
%apply SMW formula in a similar way as above to solve $\tilde p$ from the equation 
%$H(\lambda) \tilde p = M p$. Then, using \eqnok{newtonstep1}, we obtain the newtion 
%iterate $\lambda^+$ by replacing $\|w\|^2$ of \eqnok{newtonstep2} by $p^T M \tilde p$. 
%For the similar reason as discussed before, this approach also needs ${\cal O}(n)$ amount 
%of arithmetic operations for large-scale LRTR subproblems.

\subsection{Handling the hard case termination}
\label{3.2}  

Recall that one of the key parts in the implementation of steps 3) and 6) of Algorithm 1 
is the computation of a vector $u=u_\lambda$ such that $\|u\|_M=1$ and $u^TH(\lambda)u$ 
approaches zero as $\lambda \downarrow -\lambda_1$ (see Subsections \ref{2.2} and \ref{2.4}). 
In this subsection, we provide an efficient approach to find such a vector $u$ in the 
context of the low-rank version of \eqnok{subproblem}, which completely avoids
the computation of the Cholesky factorization of $H(\lambda)$.

Recall from \eqnok{hessian} that $H = D + V E V^T$, where $D \succ 0$ and $E$ 
is diagonal and nonsingular. We can partition $E$ (after performing a symmetric 
permutation of its rows and columns) as $E=\Diag(E_1, -E_2)$, where both 
$E_1$ and $E_2$ are positive diagonal matrices. Accordingly,
we partition $V$ as $V=(V_1, V_2)$, and hence
\beq \label{partition1}
V E V^T = V_1 E_1 V_1^T -  V_2 E_2 V_2^T.
\eeq
Noting that $H(\lambda) = H + \lambda M$, we can write $H(\lambda)$ as 
\beq \label{partition2}
H(\lambda) = F(\lambda) - V_2 E_2 V_2^T,
\eeq
where $F(\lambda) = D + \lambda M +  V_1 E_1 V_1^T \succ 0$ for any $\lambda \geq 0$ 
due to the fact that $D,\, M \succ 0$.

\mgap

The following technical lemma provides the key tool for our analysis in this 
subsection.

\begin{lemma} \label{hardcase}
Assume that $\lambda_1$ defined in \eqnok{lambda1} is nonpositive. Then, 
\[
\lim_{\lambda \downarrow -\lambda_1} \lambda_{\max}\left(V_2^T H(\lambda)^{-1} V_2\right) =
\infty,
\]
where $V_2$ is defined as above.
\end{lemma}

\begin{proof}
For any $\lambda > -\lambda_1$, using \eqnok{partition2} and SMW formula twice, 
we have
\beq
{H(\lambda)}^{-1} 
= {F(\lambda)}^{-1} + {F(\lambda)}^{-1} V_2 {\left(E_2^{-1} - 
{V_2}^T {F(\lambda)}^{-1} V_2\right)}^{-1} {V_2}^T {F(\lambda)}^{-1},
\label{inv1} 
\eeq
and 
\beqa
{\left(E_2^{-1} - {V_2}^T {F(\lambda)}^{-1} V_2\right)}^{-1}
&=& E_2 + E_2 {V_2}^T{\left(F(\lambda)- V_2 E_2 {V_2}^T\right)}^{-1} V_2 E_2, \nn \\
&=& E_2 + E_2 {V_2}^T {H(\lambda)}^{-1} V_2 E_2. \label{inv2} 
\eeqa 
Using the definition of $F(\lambda)$ and the fact $M \succ 0$, we easily 
see that $F (\lambda) \succ D \succ 0$ for any $\lambda > -\lambda_1 \geq 0$.
This implies that $\|F(\lambda)^{-1}\|$, and hence $\|F(\lambda)^{-1} V_2\|$, 
is bounded for all $\lambda > -\lambda_1$. From \eqnok{inv1}, we obtain that
\beq \label{normineq1}
\|{H(\lambda)}^{-1}\|
\leq \|{F(\lambda)}^{-1}\| + \|{F(\lambda)}^{-1}V_2\| \ \left\|{\left(E_2^{-1} -
    {V_2}^T {F(\lambda)}^{-1} V_2\right)}^{-1}\right\| \ \|{V_2}^T {F(\lambda)}^{-1}\|
\eeq
By the definitions of $H(\lambda)$ and $\lambda_1$, we have 
$\lim_{\lambda \downarrow -\lambda_1} \|{H(\lambda)}^{-1}\| = \infty$,
which, together with \eqnok{normineq1} and the fact that $\|F(\lambda)^{-1}\|$ 
and $\|F(\lambda)^{-1}V_2\|$ are bounded for all $\lambda > -\lambda_1 \geq 0$, 
implies  
\beq \label{limit2}
\lim_{\lambda \downarrow -\lambda_1} \left\|{\left(E_2^{-1} - {V_2}^T 
{F(\lambda)}^{-1}V_2\right)}^{-1}\right\| =\infty.
\eeq
Moreover, from \eqnok{inv2}, we have 
\beqas
\left\|{\left(E_2^{-1} - {V_2}^T {F(\lambda)}^{-1} V_2 \right)}^{-1} \right\| 
\leq \|E_2\| + \|E_2\| \ \left\|{V_2}^T {H(\lambda)}^{-1} V_2\right\| \ \|E_2\|, 
\eeqas
with together with \eqnok{limit2} implies that
\beqas 
\lambda_{\max}\left({V_2}^T {H(\lambda)}^{-1} V_2\right)= \left\|{V_2}^T 
{H(\lambda)}^{-1} V_2\right\| \to \infty \ \hbox{as} \ \lambda
\downarrow -\lambda_1. 
\eeqas \qed
\end{proof}

\vgap

The following theorem provides an efficient approach to compute the vector 
$u$ for dealing with the hard case termination in step 3) of Algorithm 1 
and updating $\lambdaL$ in step 6) of Algorithm 1.

\begin{theorem} \label{hardcase1}
Assume that $\lambda_1$ defined in \eqnok{lambda1} is nonpositive. 
Suppose $u_{\lambda} = {H(\lambda)}^{-1}v/\|{H(\lambda)}^{-1}v\|_M$,  
where $v = V_2 r$ and $r$ is a unit eigenvector of $V_2^T {H(\lambda)}^{-1} V_2$ 
corresponding to its maximum eigenvalue. 
%$\lambda_{\max}\left(V_2^T {H(\lambda)}^{-1} V_2\right)$.
Then, 
\[
\lim_{\lambda \downarrow -\lambda_1} u_{\lambda}^T H(\lambda) u_{\lambda} = 0. 
\]
\end{theorem}

\begin{proof}
%Let $v = V_2 r$. Then $u = \frac{{H(\lambda)}^{-1}v}{\|{H(\lambda)}^{-1}v\|_M}$. 
It follows from Cauchy-Schwarz inequality that 
\beq \label{cauchy}
v^T{H(\lambda)}^{-1}v \leq \|v\|_{M^{-1}} \|{H(\lambda)}^{-1}v\|_M. 
\eeq
Using \eqnok{cauchy} and the definitions of $v$ and $r$, we have
\beqas
u_{\lambda}^TH(\lambda)u_{\lambda} 
&= &  \frac{v^T H(\lambda)^{-1}v} {\|{H(\lambda)}^{-1}v\|^2_M} 
=  \frac{(v^T{H(\lambda)}^{-1}v)^2} {(v^T{H(\lambda)}^{-1}v)
\|H(\lambda)^{-1}v\|^2_M}, \\ [4.0pt]
&\leq &  \frac{\|v\|^2_{M^{-1}}}{v^T H(\lambda)^{-1} v}
= \frac{r^T V_2^T M^{-1}V_2 r}{r^T V_2^T {H(\lambda)}^{-1} V_2 r} 
, \\ [4.0pt] 
& \leq & \frac{\|V_2^T M^{-1} V_2\|}{\lambda_{\max}(V_2^T 
{H(\lambda)}^{-1} V_2)}, 
\eeqas
which, together with  Lemma \ref{hardcase}, immediately implies that 
the conclusion holds. \qed
\end{proof}      

%\vgap

Before ending this subsection, we make two observations. 
First, since $\lambda^*=-\lambda_1$ in the hard case, Theorem \ref{hardcase1} 
implies that the vector $u_{\lambda}$ defined in its statement satisfies
$\lim_{\lambda \downarrow \lambda^*}u^T_{\lambda} H(\lambda) u_{\lambda} = 0$,
which is exactly the condition required in the discussion of the
hard case (see Subsection \ref{2.2}). Second, since the number of columns of $V$ 
is assumed to be small in the low-rank version subproblem \eqnok{subproblem}, 
it follows that the matrix $V_2^T {H(\lambda)}^{-1} V_2$ is small-sized, and 
hence a unit eigenvector $r$ as in Theorem \ref{hardcase1} can be easily computed. 
Moreover, the SMW formula can be used to compute $V_2^T {H(\lambda)}^{-1} V_2$ and 
$u_{\lambda}$ in ${\cal O}(n)$ arithmetic operations. 
%Also, since $M$ has the lower-rank 
%structure as specified in \eqnok{scale}, we easily see from \eqnok{alpha} 
%that the amount of arithmetic operations for computing such a $\alpha$ in step 3  
%is ${\cal O}(n)$.

\subsection{Improving $\lambdaL$ when $\lambda < -\lambda_1$}
\label{3.3}

Recall that the key part in the implementation 
of step 4) of Algorithm 1 consists of finding a vector $u$ satisfying
$\|u\|_M=1$ and $u^TH(\lambda)u <0$, whenever $\lambda < -\lambda_1$
(see Subsections \ref{2.4}).
In this subsection, we provide an efficient approach 
to find such a vector $u$ in the context of the low-rank version
of \eqnok{subproblem}, which completely avoids the computation of a
partial Cholesky factorization of $H(\lambda)$.

Assume then that $0 \le \lambda < -\lambda_1$. This implies that
$H(\lambda)$ is indefinite, and hence that the matrix $X$
defined in \eqnok{matX1} is also indefinite, in view of
Theorem \ref{testpsd}(iii). Hence, letting
$y$ be an eigenvector of $X$ corresponding to its minimum eigenvalue,
we have that $y^TXy < 0$.
Using the definition of $\hW$ and relations \eqnok{Hlambda},
\eqnok{hDVE} and \eqnok{matX1}, we easily see that
\beq \label{reduction}
X = \hW + \hW \hE \hW = \hV^T \hD^{-1} H(\lambda) \hD^{-1} \hV.
\eeq
Hence, letting $u := \hD^{-1} \hV y/\|\hD^{-1} \hV y\|_M$, we have that
$\|u\|_M=1$ and
\[
u^T H(\lambda) u  = \frac{y^T X y}{\|\hD^{-1} \hV y\|_M^2} < 0.
\]

Note that since $X$ is a small-sized matrix (see Subsection \ref{3.1}),
it is relatively cheap to compute the vector $y$ as described above.
Moreover, since the amount of arithmetic operations to compute $X$
for a large-scale LRTR subproblem is ${\cal O}(n)$, 
the computation of the vector $u$ described above can be carried out
in ${\cal O}(n)$ arithmetic operations.

\subsection{Finding the initial $\lambdaU$}
\label{3.4}

Recall that step 1) of Algorithm 1 requires initial estimates of
the lower bound $\lambdaL$ and the upper bound $\lambdaU$.
An approach for estimating these bounds for a general TR subproblem
in ${\cal O}(n^2)$ arithmetic operations is described in Section 7.3.8
of \cite{CONN00}. For the large-scale lower-rank version subproblem
\eqnok{subproblem}, the above approach is expensive, and hence
not suitable. 

In our implementation of Algorithm 1, we set $\lambdaL = 0$.
We now provide an efficient approach to find an initial estimate of $\lambdaU$ 
in the context of the low-rank version of \eqnok{subproblem} in this 
subsection. 

Recall that $H$ and $M$ have low-rank structure (see \eqnok{hessian} and \eqnok{scale}).
Assume that the row size of matrices $E$ and $\tE$ is $\bar k$ and $\tilde k$, 
respectively. For the convenience of the presentation, we rewrite $H$ and $M$ 
in \eqnok{hessian} and \eqnok{scale} as follows
\beqa
H = D + \sum_{i=1}^{\bar k} E_{ii} v_i v_i^T, \label{hessian1} \\
M = \tD + \sum_{i=1}^{\tilde k} \tE_{ii} \tv_i \tv_i^T, \label{scale1}
\eeqa
where $v_i, \, \tv_i$ are the $i$th column of $V$ and $\tV$, respectively.

Given any $\epsilon>0$, we can trivially set initial $\lambdaU$ 
to be $\epsilon$ if $\lambda^*=0$. Hence, we now assume that $\lambda^* >0$. 
It follows from Lemma \ref{optcond} that $\lambda^*$ together with a global solution 
$p$ of \eqnok{subproblem} satisfies  
\beqa
(H + \lambda^* M) p & = & -g,  \label{bound1} \\
\|p\|_M & = &  \Delta, \label{bound2} \\
H + \lambda^* M & \succeq 0 &. \label{bound3}
\eeqa
Multiplying \eqnok{bound1} by $p^T$ on the left and using  
\eqnok{bound2} and the fact $M \succ 0$, we obtain
\[
p^T H p + \lambda^* \Delta^2 = -p^T g \leq \|p\|_M \|g\|_{M^{-1}}
= \Delta \|g\|_{M^{-1}}.
\]
Hence
\beq \label{ineq1}
\lambda^* \leq \|g\|_{M^{-1}} {\Delta}^{-1} - (p^T H p){\Delta}^{-2}.
\eeq 
Let $\tp = \Mh p / \Delta$. Noting that $\|\tp\| = 1$, we have
\beqa
(p^T H p){\Delta}^{-2}
&  =  & \tp^T\Mhn H \Mhn\tp, \nn \\
& \geq & \tp^T \Mhn \left(D + \sum\limits_{\{i|E_{ii} < 0\}}
{E_{ii} v_i v_i^T}\right)\Mhn \tp, \nn \\ [4pt]
& \geq & \lambda_{\min}(D) \|\Mhn \tp\|^2 + \sum\limits_{\{i|E_{ii} < 0\}} 
E_{ii} (v_i^T\Mhn \tp)^2, \nn \\ 
& \geq & \zeta \|\Mhn \tp\|^2, \label{lowerbdd01}
\eeqa
where
$\zeta = \lambda_{\min}(D) + \sum\limits_{\{i|E_{ii} < 0\}}{E_{ii}\|v_i\|^2}$.
Using \eqnok{scale1} and the fact that $\tE_{ii}>0$ for all $i$, we have
\beq \label{eigenM}
\lambda_{\min}(M) \ge \min\limits_{1 \le i \le n} \tD_{ii} \,\,\, \mbox{and} \,\,\,
\lambda_{\max}(M) \le \max\limits_{1 \le i \le n} \tD_{ii} + \sum\limits_{i=1}^{\tilde k}   
\tE_{ii} \|\tv_i\|^2.
\eeq
Note that
\beqas
{\lambda_{\max}(M)}^{-1} \leq \|\Mhn \tp\|^2 \leq
{\lambda_{\min}(M)}^{-1}.
\eeqas
This together with \eqnok{eigenM} implies that 
\beq \label{Mhp}
{\left(\max\limits_{1 \le i \le n} \tD_{ii} + \sum\limits_{i=1}^{\tilde k}\tE_{ii} \|\tv_i\|^2\right)}^{-1} 
\leq \|\Mhn \tp\|^2 \leq {\left(\min\limits_{1 \le i \le n} \tD_{ii}\right)}^{-1}.
\eeq 
Using \eqnok{ineq1}, \eqnok{lowerbdd01}, and \eqnok{Mhp}, we see that 
$\lambda^*  \leq \lambda^e$, where $\lambda^e$ is defined as   
\beq \label{lambdaU}
\lambda^e \equiv \left\{ \begin{array}{ll}
\|g\|_{M^{-1}}{\Delta}^{-1} -  \zeta{\left(\max\limits_{1 \le i \le n} \tD_{ii}
+ \sum\limits_{i=1}^{\tilde k} \tE_{ii} \|\tv_i\|^2 \right)}^{-1} & \mbox{if $\zeta \ge 0$},
\\ \\
\|g\|_{M^{-1}}{\Delta}^{-1} - \zeta {\left(\min\limits_{1 \le i \le n} \tD_{ii}\right)}^{-1}  &
\mbox{if $\zeta < 0$}. \end{array}
\right.
\eeq
This together with the fact that $\lambdaU=\epsilon$ if $\lambda^*=0$ implies that 
$\max(\epsilon, \lambda^e)$ is a proper initial estimate of the upper bound $\lambda^U$ 
for any $\lambda^*$.

Using the fact that $M$ \eqnok{scale} has low-rank structure, we see that $M^{-1}g$ 
can be computed by means of SMW formula requiring ${\cal O} (n)$ arithmetic operations.  
This together with \eqnok{lambdaU}, and the fact that 
$\|g\|_{M^{-1}} = \sqrt{g^T M^{-1} g}$ and $\bar k$ and $\tilde k$ are small, 
implies that this approach for finding initial estimate of the upper bound $\lambdaU$ 
requires ${\cal O} (n)$ arithmetic operations in the context of a LRTR subproblem.      
   
\section{Some numerical implementation results}
\label{4}

In this section, our main goal is to test the numerical performance of
``low-rank'' trust region methods whose LRTR
subproblems are solved by the MNE method proposed in Section \ref{3}. 
For this purpose, we implement a specific version of the modified 
log-barrier (MLB) algorithm due to Polyak \cite{POLYAK92} 
(see Subsection \ref{4.1}) and use it to solve a collection
of nonlinear programming problems from CUTEr \cite{GoOrTo03} where 
only simple bound constraints are present. Like the log-barrier 
method discussed in Section \ref{introduction}, the MLB algorithm also 
consists of solving a parametrized family of unconstrained nonlinear 
problems. In our implementation, these subproblems are solved by using 
a ``low-rank'' trust region approach similar to the one discussed in 
Section \ref{introduction} in the context of the log-barrier method.
This section is divided into two subsections. In Subsection \ref{4.1},
we discuss the generic MLB algorithm for solving nonlinear
programming problems with general inequality constraints and its
specialization to problems with simple bound constraints. 
%In Subsection \ref{4.2}, we report our computational results obtained 
%by comparing our method with LANCELOT \cite{CONN92} on the forementioned 
%collection of problems from CUTEr.
In Subsection \ref{4.2}, we report the computational results of 
our implementation of the MLB method and its comparison with a version of
LANCELOT \cite{CONN92} based on the forementioned collection of problems 
from CUTEr \cite{GoOrTo03}.

\subsection{The modified log-barrier algorithm}
\label{4.1}

In this subsection, we discuss the generic MLB algorithm for solving 
nonlinear programming problems with general inequality constraints and its
specialization to problems with simple bound constraints.
 
Consider the nonlinear programming problem
\beq \label{NLP}
\begin{array}{ll}
\mbox{minimize} & f(x) \\ [6pt]
\mbox{subject to} & c_i(x) \ge 0, \ i=1, \cdots, m.
\end{array}
\eeq
where the functions $f(x)$ and $c_i(x), \, i=1, \cdots, m$ are twice 
continuously differentiable in $\Re^n$. Its associated first-order 
optimality conditions are: 
\beqas
\nabla f(x) - \sum\limits_{i=1}^m \lambda_i \nabla c_i(x)=0, \ \lambda \ge 0, \\
\sum_{i=1}^m \lambda_i c_i(x)=0, \ c_i(x) \ge 0, \ i=1, \cdots, m.
\eeqas
The MLB method proposed by Polyak \cite{POLYAK92} 
consists of solving a sequence of unconstrained problems with objective 
functions given by   
\beq \label{augbarr}
\cM(x, \muk, \lambdak) = f(x) - \muk \sum_{i=1}^m \lambda^{(k)}_i 
\log\left(\frac{c_i(x)}{\muk}+1\right), 
\eeq
where $\lambdak \in \Re^m$ is an estimate of a Lagrange multiplier 
at a solution of \eqnok{NLP} and $\muk > 0$ is a log-barrier parameter. 
Letting $\xk$ denote a stationary point of $\cM(x, \muk, \lambdak)$, 
Polyak \cite{POLYAK92} has shown under reasonable conditions that 
there exists a threshold value $\bar{\mu}>0$ such that for any fixed 
$\mu \in (0, \bar{\mu})$, the MLB method which updates $\{\lambdak\}$ 
according to
\beq \label{chglambda}
\lambda^{(k+1)}_i = \lambda^{(k)}_i/(c_i(\xk)/\mu+1),
\eeq
generates a sequence of iterates $\{(\xk, \lambdak)\}$ which converges 
to a point satisfying the first-order optimality condition of \eqnok{NLP} 
as $k \to \infty$. 

In order to give the detailed description of the MLB method, we introduce 
the following parameters and definitions (see Breitfeld and Shanno 
\cite{BreitShanno94} and \cite{BreitShanno96}).

Let $T_1 = 10^{-4}$, $T_2 = 10^{-6}$, $\epsilon_0 = 10^{-5}$, $\sigma=0.5$ 
or $0.1$. Also, let
\beqa
r &=& -0.5 \log_{10}(T_1), \ \ \ \epsilon_k = \max(\epsilon_0, 10^{-(k+r-1)}), \nn \\ [4.0pt]
\nu^{(k)}_1 &=& \max\left\{\frac{\|\nabla\cM\|}{1+\|\xk\|}, -\min_{i=1, \cdots, m} c_i(\xk) \right\},
\nn \\ [4.0pt]
\nu^{(k)}_2 &=& \max\left\{ \nu^{(k)}_1, \frac{\sum_{i=1}^m \lambda^{(k+1)}_i|c_i(\xk)|}{1+\|\xk\|},
\frac{\|\nabla f(\xk) - \sum_{i=1}^m \lambda^{(k+1)}_i \nabla c_i(\xk)\|_{\infty}}{1+\|\xk\|} 
\right\}, \label{nk1} \\ [4.0pt]
\nu^{(k)}_3 &=& \frac{|f(\xk)-f(x^{(k-1)})|}{1+|f(x^{(k-1)})|},\ \ \
\nu^{(k)}_4 = \frac{|\nu^{(k)}_2 - \nu^{(k-1)}_2|}{\nu^{(k-1)}_2}. \nn
\eeqa

The algorithm below is a complete description of MLB method as implemented 
in \cite{BreitShanno94} and \cite{BreitShanno96}.

\begin{tabbing}  
\={\bf Algorithm MLB:} \\ [4.0pt]
\> Le\=t $\xi \in (0,1)$, $\mu^{(1)}>0$, $\lambda^{(1)}>0$, $\nu^{(0)}_1$, $\nu^{(0)}_2$,
$\nu^{(0)}_3$ be given, and set $k = 1$. \\
\> {\bf For} \= $k = 1, 2, 3, \ldots$, until $\nu^{(k-1)}_2 < T_1$ or 
($\nu^{(k-1)}_1 < T_1$ and $\nu^{(k-1)}_3 < T_2$),
{\bf do}
\\
\> \> 1. \= Minimize the function \eqnok{augbarr} approximately, \\ 
\> \>    \> obtaining $\xk$ such that $\|\nabla\cM(\xk,\muk,\lambdak)\| \le \epsilon_k$. \\
\> \> 2. \= Update $\lambda^{(k+1)}$ by \eqnok{chglambda}. If $\nu^{(k)}_4 < \xi$, 
	    set $\mu^{(k+1)} = \sigma \mu^{(k)}$. \\
\> \>    \> Increment $k$ by 1, and return to step 1. \\
\> {\bf End}
\end{tabbing}

In our implementation, we initialize $\mu^{(0)}=10^{-2}$ and $\lambda^{(0)}_i=1$ for
$i=1, \cdots, m$ as suggested in \cite{BreitShanno94} and \cite{POLYAK92}, respectively.

While Breitfeld and Shanno \cite{BreitShanno94} 
have used a line search method to find $\xk$, we instead use a trust region method 
whose associated TR subproblems have quadratic objective functions 
with low-rank Hessian matrices obtained by means of limited-memory BFGS (L-BFGS) method 
(see Section 9.1 of \cite{NOCEDAL99})). 
The MNE method is then used to solve the resulting LRTR subproblems. 

We next provide more details of how a low-rank approximation of
the Hessian of $\cM(\cdot,\muk,\lambdak)$ is computed 
in the context of solving nonlinear programming problems with
simple bound constraints. For the purpose of this discussion,
assume that the constraints of \eqnok{NLP} are given by
$c^l_i(x) \equiv x_i - l_i \geq 0$ for $i \in I_l$ and
$c^u_i(x) \equiv u_i- x_i \geq 0$ for $i \in I_u$
where $I_l, I_u \subseteq \{1,\cdots,n\}$ are index sets corresponding
to the lower and upper bound constraints, respectively.
%If a problem is unconstrained, we set $I_l = I_u = \emptyset$.
Let $\lambda^l$ and $\lambda^u$ be the Lagrange multipliers
corresponding to the lower and upper bound constraints, respectively. 
We easily see from \eqnok{augbarr} that 
$\nabla^2 \cM(x,\mu,\lambda)=\nabla^2 f(x) + \mu Q$, where
\beq \label{aughess}
Q = \sum_{i \in I_l} \frac{\lambda^l_ie_i e^T_i}{(c^l_i(x)+\mu)^2}
+ \sum_{i \in I_u} \frac{\lambda^u_ie_i e^T_i}{(c^u_i(x)+\mu)^2}.
\eeq
Since $Q$ is a diagonal matrix and can be easily computed, it makes
sense to just compute a low-rank approximation $F$ of
$\nabla^2 f(x)$ and use $H \equiv F + \mu Q$ as a low-rank
approximation of $\nabla^2 \cM(x,\mu,\lambda)$.
We use the L-BFGS method to obtain a low-rank apprximation of
$\nabla^2 f(x)$. This matrix $H$ is used as the Hessian of the
objective function of the LRTR subproblem \eqnok{subproblem}
at the point $x$.  Moreover, the other data of this subproblem
is determined as $g \equiv \nabla \cM (x, \mu, \lambda)$ and
$M \equiv \Diag(v)$, where $v \in \Re^n$ is defined as
$v_i= 1/\sqrt{Q_{ii}}$ if $i \in I_l \cup I_u$, and $v_i=1$
if $i \notin I_l \cup I_u$.

\subsection{Implementations of some problems from CUTEr}
\label{4.2}

%In this subsection, we will report the computational performance of our implementation 
%of a version of Algorithm MLB described in Subsection \ref{4.1} for solving a collection 
%of problems from CUTEr \cite{GoOrTo03} with no constraint or some simple bound constraints. 
%We also compare the computational results obtained by our method with those obtained by 
%LANCELOT \cite{CONN92} for those test problems.  
In this subsection, we will report the computational results obtained from the
implementation of a specifc version of Algorithm MLB described in Subsection \ref{4.1} 
and present the comparisons of our method with LANCELOT \cite{CONN92} on a collection
of nonlinear programming problems from CUTEr \cite{GoOrTo03} where
only simple bound constraints are present.

All computations are performed on a Sun Ultra 10 workstation which has a
single UltraSPARC IIi processor running at 440Mhz and 512MB of memory.
The sixty test problems are selected from CUTEr. Seventeen of them have simple 
bound constraints and fixed variables. The remaining problems are unconstrained ones. 
Each row of Tables \ref{problems-cuter1} and \ref{problems-cuter2} gives the problem name, 
the number of variables, the number of bound constraints, 
the number of free variables, and the number of fixed variables
on columns one through five, respectively. 

%In the discussion below, we will refer to  
Some computational results for our code are also presented in 
Tables \ref{problems-cuter1} and \ref{problems-cuter2}. For each test problem, 
the total number of iterations performed by the MNE method is given in the sixth column, 
and the total number of LRTR subproblems \eqnok{subproblem} solved is given in the seventh 
column. The average iterations performed by the MNE method for each test problem 
is given in the eighth column, which is obtained by dividing the entry in the sixth column
by the entry in the seventh clolumn.   

 From the eighth column of Tables \ref{problems-cuter1} and \ref{problems-cuter2}, 
we observe that the average iterations performed by the MNE method are between 1.0 and 2.0 
for almost all test problems, except the problems QRTQUAD (2.89) and FLETCBV2 (0). 
Moreover, the average of the entries in the eighth column over all test problems is $1.25$, 
which is better than the average $1.6$ obtained by Mor{\'e} and Sorensen \cite{MORE83} 
for the NE method applied to the general TR subproblems \eqnok{subproblem}. Hence, 
this indicates that the MNE method developed in our paper is efficient and robust for 
solving the low-rank version TR subproblem \eqnok{subproblem}.  

\begin{table}[t]
\caption{The main test problems and some results of MTR}
\centering
\label{problems-cuter1}
%\begin{center}
\begin{small}
\begin{tabular}{lrrrrrrr}
\multicolumn{1}{c}{Problem} & \multicolumn{1}{c}{Var} &
\multicolumn{1}{c}{Free} & \multicolumn{1}{c}{Bound} &
\multicolumn{1}{c}{Fix} & \multicolumn{1}{c}{TrsIt} & 
\multicolumn{1}{c}{Trs} & \multicolumn{1}{c}{Ratio} \\
\hline
ARGLINA   &    387  &    387  &     0  &   0  & 11     & 6       & 1.83\\ 
ARGLINB   &    387  &    387  &     0  &   0  & 25     & 21      & 1.19\\
ARGLINC   &    388  &    388  &     0  &   0  & 23     & 19      & 1.21 \\
ARWHEAD   &  10000  &  10000  &     0  &   0  & 13     & 12      & 1.08 \\
BDQRTIC   &  19999  &  19999  &     0  &   0  & 180    & 136     & 1.32\\
BRATU1D   &  13333  &  13331  &     0  &   2  & 6508   & 6021    & 1.08 \\
BROWNAL   &    547  &    547  &     0  &   0  & 11     & 8       & 1.38\\
BRYBND    &  14288  &  14288  &     0  &   0  & 70     & 61      & 1.15\\ 
CHEBYQAD  &    316  &      0  &   316  &   0  & 383    & 342     & 1.12 \\
COSINE    &  20000  &  20000  &     0  &   0  & 27     & 14      & 1.93 \\
CURLY20   &  14295  &  14295  &     0  &   0  & 1597   & 1471    & 1.09 \\
CVXBQP1   &  20000  &      0  & 20000  &   0  & 1200   & 1162    & 1.03\\
DIXMAANA  &  19998  &  19998  &     0  &   0  & 21     & 15      & 1.4\\
DIXMAANB  &  19998  &  19998  &     0  &   0  & 21     & 15      & 1.4\\   
DIXMAANC  &  19998  &  19998  &     0  &   0  & 23     & 17      & 1.35\\
DIXMAAND  &  19998  &  19998  &     0  &   0  & 25     & 19      & 1.32\\
DIXMAANJ  &  19998  &  19998  &     0  &   0  & 626    & 617     & 1.01\\
DIXMAANK  &  19998  &  19998  &     0  &   0  & 419    & 413     & 1.01\\
DIXMAANL  &  19998  &  19998  &     0  &   0  & 450    & 444     & 1.01\\
EDENSCH   &  13333  &  13333  &     0  &   0  & 55     & 43      & 1.28\\  
EIGENBLS  &   1122  &   1122  &     0  &   0  & 1695   & 1608    & 1.05\\
EIGENCLS  &   1122  &   1122  &     0  &   0  & 1559   & 1452    & 1.07\\  
ENGVAL1   &  19999  &  19999  &     0  &   0  & 32     & 26      & 1.23\\  
ERRINROS  &     50  &     50  &     0  &   0  & 819    & 772     & 1.06\\
FLETCBV2  &  13330  &  13330  &     0  &   0  & 0      & 0       & 0.00 \\
FMINSRF2  &  19881  &  19881  &     0  &   0  & 846    & 841     & 1.01\\
FREUROTH  &  19999  &  19999  &     0  &   0  & 97     & 61      & 1.59\\
FMINSURF  &  19881  &  19881  &     0  &   0  & 2405   & 2362    & 1.02\\
GRIDGENA  &   6218  &      0  &  5560  & 658  & 111    & 98      & 1.13\\
HILBERTB  &    282  &    282  &     0  &   0  & 16     & 11      & 1.45\\
LIARWHD   &  19999  &  19999  &     0  &   0  & 43     & 32      & 1.34\\
LINVERSE  &  10001  &   5000  &  5001  &   0  & 110    & 85      & 1.29\\
LMINSURF  &  19881  &  19321  &     0  & 560  & 1672   & 1669    & 1.00\\
MANCINO   &    183  &    183  &     0  &   0  & 32     & 21      & 1.52\\
NLMSURF   &  19881  &  19321  &     0  & 560  & 3618   & 3579    & 1.01\\
NOBNDTOR  &  12544  &   6050  &  6050  & 444  & 413    & 346     & 1.19\\
NONDIA    &  16666  &  16666  &     0  &   0  & 19     & 18      & 1.06\\
\hline
\end{tabular}
\end{small}
%\end{center}
\end{table}

\begin{table}[t]
\caption{The main test problems and some results of MTR(cont'd)}
\centering
\label{problems-cuter2}
%\begin{center}
\begin{small}
\begin{tabular}{lrrrrrrr}
\multicolumn{1}{c}{Problem} & \multicolumn{1}{c}{Var} &
\multicolumn{1}{c}{Free} & \multicolumn{1}{c}{Bound} &
\multicolumn{1}{c}{Fix} & \multicolumn{1}{c}{TrsIt} &
\multicolumn{1}{c}{Trs} & \multicolumn{1}{c}{Ratio} \\
\hline
OBSTCLAE  &  12769  &      0  & 12321  & 448  & 296     & 254     & 1.17\\
OBSTCLBL  &  12769  &      0  & 12321  & 448  & 135     & 113     & 1.19\\
OBSTCLBM  &  12769  &      0  & 12321  & 448  & 139     & 120     & 1.16\\
OBSTCLBU  &  12769  &      0  & 12321  & 448  & 141     & 117     & 1.21\\
ODC       &  16900  &  16384  &     0  & 516  & 682     & 682     & 1.00\\                
PENALTY1  &  19999  &  19999  &     0  &   0  & 151     & 116     & 1.30\\
PENALTY3  &    120  &    120  &     0  &   0  & 94      & 76      & 1.24\\
POWELLSG  &  20000  &  20000  &     0  &   0  & 39      & 31      & 1.26\\
POWER     &  20000  &  20000  &     0  &   0  & 779     & 740     & 1.05\\
PROBPENL  &  19999  &      0  & 19999  &   0  & 10      & 6       & 1.67\\
QRTQUAD   &  20000  &  10000  & 10000  &   0  & 104     & 36      & 2.89\\
%RAYBENDL  &  20000  &  19996  &     0  &   4 \\
%RAYBENDS  &  20000  &  19996  &     0  &   4  \\
SENSORS   &    199  &    199  &     0  &   0  & 51      & 42      & 1.21\\    
SINQUAD   &  19999  &  19999  &     0  &   0  & 82      & 54      & 1.52\\ 
SPARSQUR  &  16666  &  16666  &     0  &   0  & 41      & 37      & 1.11\\
TOINTGSS  &  20000  &  20000  &     0  &   0  & 29      & 21      & 1.38\\
TORSION2  &  12544  &      0  & 12100  & 444  & 357     & 298     & 1.20 \\
TORSION4  &  12544  &      0  & 12100  & 444  & 203     & 180     & 1.13\\   
TORSION6  &  12544  &      0  & 12100  & 444  & 149     & 125     & 1.19\\
TORSIONB  &  12544  &      0  & 12100  & 444  & 290     & 269     & 1.08\\
TORSIOND  &  12544  &      0  & 12100  & 444  & 244     & 200     & 1.22\\
TORSIONF  &  12544  &      0  & 12100  & 444  & 143     & 123     & 1.16\\
TQUARTIC  &  19999  &  19999  &     0  &   0  & 29      & 21      & 1.38\\
VARDIM    &  19999  &  19999  &     0  &   0  & 101     & 96      & 1.05\\
\hline
\end{tabular}
\end{small}
%\end{center}
\end{table}

%All computations are performed on a Sun Ultra 10 workstation which has a 
%single UltraSPARC IIi processor running at 440Mhz and 512MB of memory.  
Our code MTR is written in ANSI C and LANCELOT is a FORTAN 
code. MTR and LANCELOT are both compiled under the default optimization.  
In MTR, we set the parameter $\sigma = 0.1$ for the Algorithm MLB and store
3 most recent vector pairs that provide curvature information for L-BFGS update. 
We select the same initial point $x^0$ as LANCELOT 
whenever it is strictly feasible; otherwise, we modify the infeasible components of 
this initial point $x^0$ in order to make it strictly feasible. We set up an upper 
bound of one hour computation 
time (or 3,600 seconds) per problem for both codes. We implement the version of LANCELOT 
which uses bfgs-approximate-second-derivatives, bandsolver-preconditioned-cg-solver, 
inexact-cauchy-point, two-norm-trust-region and all its other default settings. 
LANCELOT terminates when the infinity norm of the projected gradient is less than 
$10^{-5}$ or exceeds one hour computation time. 

Tables \ref{compres-cuter1} and 
\ref{compres-cuter2} give the performance of MTR and LANCELOT. The objective function 
values of both methods are given in the second and third columns.
The CPU times (in seconds) are given in the fourth and fifth columns.
The iterations given in the sixth and
seventh columns represent the total number of conjugate gradient        
iterations performed by LANCELOT and the total number of TR
subproblems generated by MTR, respectively. 
The total numbers of function and gradient evaluations for both codes are given in the 
last two columns. 

%From our implementation, we observe that the average inner iterations of 
%solving the specialized TR subproblem by using the modified nearly exact method is 
%$1.33$, which is even better than $1.6$ reported in Mor{\'e} and Sorensen 
%\cite{MORE83} by the nearly exact method. Hence, the modified nearly exact 
%method developed in our paper is efficient and robust for solving the 
%specialize TR subproblem. 

Based on the results of Tables \ref{compres-cuter1} and \ref{compres-cuter2}, we now give some 
conclusions about the performance of the MTR and LANCELOT codes 
in terms of CPU time and the relative difference (rel. diff.) of the objective values obtained by 
MTR and LANCELOT. 
Among those test problems, MTR has:
\begin{itemize}
\item[i)]  
better CPU time and better or close optimal value (i.e., rel. diff. $\leq$ 1.0e-5) for 
$60 \%$ of problems;
\item[ii)]
better CPU time and worse optimal value (i.e., rel. diff. $>$ 1.0e-5) for $10 \%$ of problems;       
\item[iii)]
worse CPU time and better or close optimal value for $23 \%$ of problems;
\item[iv)]
worse CPU time and worse optimal value for $7 \%$ of problems.
\end{itemize}
Based on the above comparison, we see that it is promising to 
solve large-scale problems with simple bound constraints
using our approach.

\begin{table}[t]
\caption{Comparison of the Two Methods on the main test problems}
\centering
\label{compres-cuter1}
%\begin{center}
\begin{small}
\begin{tabular}{lrrrrrrrr}
\multicolumn{1}{c}{Problem} & \multicolumn{2}{c}{Obj Value} &   
\multicolumn{2}{c}{Time} & \multicolumn{2}{c}{Iter} &
\multicolumn{2}{c}{Nfg} \\
\multicolumn{1}{c}{Name} & \multicolumn{1}{c}{\sc mtr} &
\multicolumn{1}{c}{\sc lan} & \multicolumn{1}{c}{\sc mtr} &
\multicolumn{1}{c}{\sc lan} & \multicolumn{1}{c}{\sc mtr} &
\multicolumn{1}{c}{\sc lan} & \multicolumn{1}{c}{\sc mtr} &
\multicolumn{1}{c}{\sc lan} \\
\hline
ARGLINA   &  3.870000e+02  &  3.870000e+02  &     0.78 &    38.48 &      6 &    3 &    14 &    14 \\  %1 
ARGLINB   &  1.931252e+02  &  1.931252e+02  &     2.25 &    29.00 &     21 &    1 &    44 &    13 \\  %1
ARGLINC   &  1.951252e+02  &  1.951252e+02  &     2.05 &    28.89 &     19 &    2 &    40 &    13 \\  %1
ARWHEAD   &  0.000000e+00  &  0.000000e+00  &     5.52 &     7.86 &     12 &    1 &    26 &    12 \\  %1
BDQRTIC   &  8.008640e+04  &  8.008640e+04  &    82.48 &    21.89 &    136 &   17 &   274 &    40 \\  %2
BRATU1D   &  5.840393e+06  &  1.261025e+08  &  3600.00 &  3600.00 &   6021 & 4540 & 12044 & 16852 \\  %1
BROWNAL   &  2.572063e-11  &  2.615568e-11  &     0.98 &    81.35 &      8 &    5 &    18 &    14 \\  %1
BRYBND    &  1.325909e-12  &  9.511770e-14  &    26.00 &    22.94 &     60 &   54 &   122 &    78 \\  %4
CHEBYQAD  &  6.778464e-03  &  8.047224e-03  &   224.84 &  3600.00 &    342 & 2394 &   698 &  2103 \\  %1   
COSINE    & -2.000000e+04  & -1.999900e+04  &     8.35 &     5.16 &     14 &   11 &    30 &    33 \\  %2
CURLY20   & -1.434017e+06  & -1.434021e+06  &   580.10 &  3600.00 &   1471 &43104 &  2944 &    47 \\  %3
CVXBQP1   &  2.894382e-04  &  9.000450e+06  &   577.99 &   238.05 &   1162 & 9760 &  2338 &    14 \\  %2
DIXMAANA  &  1.000000e+00  &  1.000000e+00  &     8.13 &     7.83 &     15 &   21 &    32 &    34 \\  %1
DIXMAANB  &  1.000000e+00  &  1.000000e+00  &     8.46 &     9.14 &     15 &   17 &    32 &    40 \\  %1      
DIXMAANC  &  1.000000e+00  &  1.000000e+00  &     9.51 &     9.79 &     17 &   20 &    36 &    42 \\  %1
DIXMAAND  &  1.000000e+00  &  1.000000e+00  &    10.57 &    14.96 &     19 &   33 &    40 &    66 \\  %1
DIXMAANJ  &  1.000001e+00  &  1.000000e+00  &   325.88 &    52.51 &    617 &  577 &  1236 &    88 \\  %2
DIXMAANK  &  1.000001e+00  &  1.000000e+00  &   217.59 &    46.26 &    413 &  389 &   828 &    87 \\  %2
DIXMAANL  &  1.000001e+00  &  1.000000e+00  &   233.88 &    38.11 &    444 &  354 &   890 &    84 \\  %2
EDENSCH   &  8.000128e+04  &  8.000128e+04  &    15.97 &    10.18 &     43 &   25 &    88 &    70 \\  %2
EIGENBLS  &  1.974164e-03  &  1.125255e-01  &   172.40 &  3600.00 &   1608 & 3424 &  3218 &  6914 \\  %1
EIGENCLS  &  2.454111e-02  &  7.176403e+02  &   158.73 &  3600.00 &   1452 & 3131 &  2906 &  9172 \\  %1
ENGVAL1   &  2.219933e+04  &  2.219933e+04  &    14.35 &     5.83 &     26 &   10 &    54 &    28 \\  %2
ERRINROS  &  3.990416e+01  &  3.990415e+01  &     0.83 &     0.17 &    772 &  110 &  1546 &   209 \\  %2
FLETCBV2  & -5.001006e-01  & -5.001006e-01  &     0.17 &     0.59 &      0 &    0 &     2 &     2 \\  %1
FMINSRF2  &  1.000002e+00  &  1.000000e+00  &   413.21 &   676.48 &    841 & 1454 &  1684 &  2341 \\  %1
FREUROTH  &  2.433123e+06  &  2.433123e+06  &    39.90 &    11.97 &     61 &   20 &   124 &    47 \\  %2
FMINSURF  &  1.000002e+00  &  7.466224e+00  &  1176.94 &  3600.00 &   2362 &  106 &  4726 &   211 \\  %1
GRIDGENA  &  2.352000e+04  &  2.352000e+04  &    16.89 &    10.70 &     98 &  126 &   210 &    48 \\  %2
HILBERTB  &  6.262020e-16  &  8.799383e-12  &     1.41 &    12.19 &     11 &   31 &    24 &    86 \\  %1
LIARWHD   &  2.130882e-13  &  5.932646e-14  &    17.03 &    13.34 &     32 &   43 &    66 &    64 \\  %4
LINVERSE  &  3.409000e+03  &  3.409000e+03  &    29.54 &  3600.00 &     85 &12067 &   184 & 16626 \\  %1
LMINSURF  &  9.000355e+00  &  9.000000e+00  &   856.11 &  3059.60 &   1669 & 7650 &  3340 & 10168 \\  %3        
MANCINO   &  2.117119e-16  &  9.787775e-20  &    17.58 &    22.76 &     21 &   12 &    44 &    45 \\  %3
NLMSURF   &  3.914874e+01  &  4.290681e+01  &  1760.28 &  3600.00 &   3579 & 7237 &  7160 & 11997 \\  %1
NOBNDTOR  & -4.418497e-01  & -4.418594e-01  &   130.18 &    32.54 &    346 &  402 &   706 &    72 \\  %4
NONDIA    &  1.344590e-15  &  3.812729e-19  &     7.93 &     2.63 &    18 &     4 &    38 &    12 \\  %4
\hline
\end{tabular}
\end{small}
%\end{center}
\end{table}

\begin{table}[t]
\caption{Comparison of the Two Methods on the main test problems(cont'd)}
\centering
\label{compres-cuter2}
%\begin{center}
\begin{small}
\begin{tabular}{lrrrrrrrr}
\multicolumn{1}{c}{Problem} & \multicolumn{2}{c}{Obj Value} &
\multicolumn{2}{c}{Time} & \multicolumn{2}{c}{Iter} &
\multicolumn{2}{c}{Nfg} \\
\multicolumn{1}{c}{Name} & \multicolumn{1}{c}{\sc mtr} &
\multicolumn{1}{c}{\sc lan} & \multicolumn{1}{c}{\sc mtr} &
\multicolumn{1}{c}{\sc lan} & \multicolumn{1}{c}{\sc mtr} &
\multicolumn{1}{c}{\sc lan} & \multicolumn{1}{c}{\sc mtr} &
\multicolumn{1}{c}{\sc lan} \\
\hline
OBSTCLAE  &  1.894773e+00  &  1.894763e+00  &   105.41 &   389.08 &   254 &  6264 &   522 &    26 \\ %1 
OBSTCLBL  &  7.285840e+00  &  7.285834e+00  &    45.48 &   165.98 &   113 &  2938 &   240 &    36 \\ %1
OBSTCLBM  &  7.285838e+00  &  7.285834e+00  &    48.09 &   291.40 &   120 &  4431 &   254 &    16 \\ %1
OBSTCLBU  &  7.285840e+00  &  7.285834e+00  &    47.26 &    58.55 &   117 &  1107 &   248 &    40 \\ %1
ODC       & -1.137898e-02  & -1.082613e-02  &   374.75 &  3600.00 &   682 &  6018 &  1366 & 14458 \\ %1
PENALTY1  &  1.985763e-01  &  1.985846e-01  &    61.19 &  1769.96 &   116 &    30 &   206 &   124 \\ %1
PENALTY3  &  9.998808e-04  &  7.704032e-02  &     5.63 &  3600.00 &    76 & 14981 &   154 & 25899 \\ %1
POWELLSG  &  2.023535e-07  &  3.498944e-05  &    13.52 &     4.90 &    31 &    19 &    64 &    40 \\ %2
POWER     &  2.723993e-07  &  1.546668e-08  &   305.59 &  1222.89 &   740 &    68 &  1482 &    80 \\ %3
PROBPENL  &  1.000000e-08  &  1.007191e-08  &     4.03 &   202.40 &     6 &     5 &    22 &    18 \\ %1   
QRTQUAD   & -5.375101e+10  & -1.526995e+10  &    34.46 &  3600.00 &    36 &  7354 &    72 &    23 \\ %1
SENSORS   & -8.336250e+03  & -8.064563e+03  &    18.57 &   265.23 &    42 &   282 &    86 &  1008 \\ %1 
SINQUAD   & -1.040172e+08  &  1.509966e-05  &    34.91 &   208.92 &    54 &   781 &   110 &   586 \\ %1
SPARSQUR  &  9.098914e-10  &  3.009349e-06  &    17.47 &    19.87 &    37 &    68 &    76 &    46 \\ %1
TOINTGSS  &  1.000050e+01  &  1.000050e+01  &    11.79 &     8.34 &    21 &    23 &    44 &    48 \\ %2  
TORSION2  & -4.263058e-01  & -4.263350e-01  &   117.85 &   340.90 &   298 &  5327 &   610 &    28 \\ %3
TORSION4  & -1.212871e+00  & -1.212881e+00  &    69.53 &   459.81 &   180 &  7990 &   374 &    20 \\ %1%2 
TORSION6  & -2.859436e+00  & -2.859446e+00  &    49.44 &   474.08 &   125 &  9381 &   264 &    18 \\ %1  
TORSIONB  & -4.183918e-01  & -4.184089e-01  &   106.21 &   267.73 &   269 &  3829 &   552 &    20 \\ %3
TORSIOND  & -1.204444e+00  & -1.204454e+00  &    82.87 &   497.59 &   200 &  8585 &   414 &    22 \\ %1%2
TORSIONF  & -2.850759e+00  & -2.850769e+00  &    50.39 &   485.76 &   123 &  9229 &   260 &    16 \\ %1%2
TQUARTIC  &  1.098962e-14  &  2.620727e-09  &    10.51 &    15.01 &    21 &    53 &    44 &    50 \\ %1
VARDIM    &  1.524712e-15  &  3.050031e-04  &    40.28 &  3600.00 &    96 &     2 &   194 &   108 \\ %1
\hline 
\end{tabular}
\end{small}
%\end{center}
\end{table}

\begin{acknowledgements} 
The authors would like to thank Nicholas Gould, Jorge Nocedal, 
Dominique Orban and Philippe Toint for their kind help in the process
of installing the softwares CUTEr and LANCELOT on our machine.
\end{acknowledgements}

\begin{thebibliography}{3}

\bibitem{BreitShanno94}
Breitfeld, M.G., Shanno, D.F.: 
Preliminary computational experience with modified log-barrier functions
for large-scale nonlinear programming. In: Hager, W., Hearn, D., Pardalos, P. 
(ed.) Large Scale Optimization: State of the Art, 45--67, Kluwer Academic 
Publishers, Dordrecht, Boston (1994)

\bibitem{BreitShanno96}
Breitfeld, M.G., Shanno, D.F.: 
Computational experience with penalty-barrier methods for nonlinear programming. 
Annals of Operations Research {\bf 62}, 439--463 (1996)

\bibitem{BunPar71}
Bunch, J.R. and Parlett, B.N.:
Direct methods for solving symmetric indefinite systems of linear
equations. SIAM Journal on Numerical Analysis {\bf 8}, 639--655 (1971)

\bibitem{ClMoStWi79}
Cline, A.K., Moler, C.B., Stewart, G.W., Wilkinson, J.H.:
An estimate for the condition number of a matrix.
SIAM Journal on Numerical Analysis {\bf 16}, 368--375 (1979) 

\bibitem{CONN92}
Conn, A.R., Gould, N.I.M., Toint, P.L.:
LANCELOT: a Fortran package for large-scale nonlinear optimization.
Springer Series in Computational Mathematics {\bf 17},
Springer Verlag, Heidelberg, New York (1992)

\bibitem{CONN00}
Conn, A.R., Gould, N.I.M., Toint, P.L.:
Trust-region methods. SIAM Publications, Philadelphia, 
Pennsylvania, USA (2000)

\bibitem{DeMe79}
Dennis, J.E., Mei, H.W.:
Two new unconstrained optimization algorithms which use function and 
gradient values. Journal of Optimization Theory and Application {\bf 28}, 
453--482 (1979)

\bibitem{FLETCHER80}
Fletcher, R.:
Practical Methods of Optimization. Unconstrained Optimization {\bf 1}, 
John Wiley, New York (1980)

\bibitem{FORTIN02}
Fortin, C., Wolkowicz, H.:
The trust region subproblem and semidefinite programming.
Optimization methods and Software {\bf 19} (1), 41--67 (2004)

\bibitem{GAY81}
Gay, D.M.:
Computing optimal locally constrained steps.
SIAM Journal on Scientific and Statistical Computing {\bf 4} (2), 
186--197 (1981)

\bibitem{GolVL}
Golub, G.H., Van Loan, C.E.:
Matrix Computations: Second Edition.
The John Hopkins University Press, Baltimore, MD 21211 (1989)

\bibitem{GoLuRoTo99}
Gould, N.I.M., Lucidi, S., Roma, M., Toint, P.L.:
Solving the trust-region subproblem using the Lanczos method.
SIAM Journal on Optimization {\bf 9} (2), 504--525 (1999)

\bibitem{GoOrTo03}
Gould, N.I.M., Orban, D., Toint, P.L.:
General CUTEr documentation. Technical Report TR/PA/02/13, CERFACS, 
Toulouse, France (2003)

\bibitem{Hebden73}
Hebden, M.D.:
An Algorithm for minimization using exact second derivatives.
Atomic Energy Research Establishment, Report T.P.515, 
Harwell, England (1973)

\bibitem{MORE83}
Mor{\'e}, J.J., Sorensen, D.C.:
Computing a trust region step. 
SIAM Journal on Scientific and Statistical Computing {\bf 4} (3), 
553--572 (1983)

\bibitem{MORE831}
Mor{\'e}, J.J.:
Recent developments in algorithms and software for trust region 
methods. In: Bachem, A., Grotschel, M., Korte B. (ed.) 
Mathematical programming: the state of the art, University of Bonn, 
258--287, Springer-Verlag, Berlin, New York (1983)

\bibitem{NOCEDAL99}
Nocedal, J., Wright, S.J.:
Numerical optimization. Springer-Verlag, New York, USA (1999)
  
\bibitem{POLYAK92}
Polyak, R.:
Modified barrier functions (theory and methods).
Mathematical Programming {\bf 54}, 177--222 (1992)

\bibitem{POWELL70} 
Powell, M.J.D.:
A hybrid method for nonlinear equations.
In: Rabinowitz, P. (ed.) 
Numerical Methods for Nonlinear Algebraic Equations, Gordon and Breach, 
New York (1970)   

\bibitem{Reinsch71}
Reinsch, C.H.:
Smoothing by spline functions II.
Numerical Mathematics {\bf 16}, 451--454 (1971) 

\bibitem{ShScBy85}
Shultz, G.A., Schnabel, R.B., Byrd, R.H.:
A family of trust-region-based algorithms for unconstrained minimization 
with strong global convergence properties.
SIAM Journal on Numerical Analysis {\bf 22}, 47--67 (1985)

\bibitem{SORENSEN82a}
Sorensen, D.C.:
Newton's method with a model trust region modification.
SIAM Journal on Numerical Analysis {\bf 19}, 404--426 (1982)

\bibitem{STEIHAUG83}
Steihaug, T.:
The conjugate gradient method and trust regions in large-scale optimization.
SIAM Journal on Numerical Analysis {\bf 20}, 626--637 (1983)

\bibitem{Toint81}
Toint, P.L.:
Towards an efficient sparsity exploiting newton method for minimization. 
In: Duff, I.S. (ed.) Sparse matrices and their uses, 
Institute of Mathematics and its Applications Conference Series, xii+387, 
London, Academic Press Inc. (1981)

\bibitem{WaWeYu04}
Wang, Z.H., Wen, Z.W., Yuan, Y.:
A subspace trust region method for large scale unconstrained optimization. 
In: Yuan, Y. (ed.) Numerical Linear Algebra and Optimization, 
Proceeding of the 2003 International Conference on Numerical Optimization and 
Linear Algebra, 265--274, Science Press, Beijing, New York (2004) 

\bibitem{Yuan04}
Yuan, Y.:
A Subspace Trust Region Algorithm. Presented at the Conference of Multiscale 
Optimization Methods and Applications, Center for Applied Optimization, 
University of Florida, USA, February 26-28, 2004.

\bibitem{ZhXu99}
Zhang, J., Xu, C.:
A class of indefinite dogleg path methods for unconstrained minimization.
SIAM Journal on Optimization {\bf 9} (3), 646--667 (1999)
  
\end{thebibliography} 
\end{document}
