
\chapter{PATH 4.x}\label{chap:path}

% This chapter is an update to sections of the following papers:
%   ``Feasible Descent Algorithms''
%   ``Interfaces to PATH 3.0''
%   ``Complementarity in GAMS and the PATH solver''
% Updated documentation for the latest version of the code

Several algorithms for solving mixed complementarity problems, such as
MILES \cite{rutherford:miles} and PATH 
\cite{dirkse:robust,dirkse.ferris:path,dirkse.ferris:pathsearch}, are based 
upon successive linearizations: at each iteration, a linear complementarity 
problem is solved to find the next iterate.  One difficulty frequently 
encountered by these methods in practice is what to do when the linear 
subproblem has no solution, a situation that the theory associated with these 
methods \cite{ralph:global} assumes does not occur.  This chapter proposes a
theoretically justifiable escape mechanism that uses a differentiable merit 
function in conjunction with the direction generated by a linearization 
method.  The resulting framework forms the basis for PATH 4.x.

Merit functions are used extensively in the development of globalization
theory and the implementation of robust algorithms.  Broadly speaking,
a merit function summarizes how close the current iterate is to a solution 
of the problem under consideration with a single number.  In complementarity 
problems and nonlinear systems of equations, the merit functions are normally 
nonnegative, and zero precisely at a solution to the original problem.
Each merit function is typically used in a globalization strategy that
involves searching between the current iterate and the Newton point (the 
solution of the linearization).   

The classical example of a merit function in nonlinear equation solving 
is the square of the two-norm residual that measures the sum of squares of the
errors in satisfying the equations.  This merit function has one additional 
property to those listed above: namely, it is everywhere differentiable 
provided that the equation itself is everywhere differentiable. 
When the linear subproblem cannot be solved, a gradient-based descent 
direction for the differentiable merit function can be searched, 
guaranteeing that progress is made toward a stationary point of the merit 
function.  

In complementarity, the two classical merit functions are based on
the natural residual \cite{luo.tseng:error} and the normal map 
\cite{robinson:normal}.  Both the natural residual and the normal map 
provide reformulations of the complementarity problem as a system of 
equations; unfortunately, the systems and corresponding residual merit 
functions are nonsmooth.  Even with this drawback, \cite{ralph:global} 
showed how to construct an extension of the line search procedure for 
smooth nonlinear equations that enables fast local convergence of 
linearization methods based on the normal map \cite{robinson:newtons} 
under conditions that are the exact generalizations of those required for 
smooth systems.  PATH 2.9 
\cite{dirkse:robust,dirkse.ferris:path,dirkse.ferris:pathsearch} implements 
this method and uses heuristics when the linear subproblem cannot be solved.  
While these heuristics are quite successful in practice, this situation is 
nonetheless unsatisfactory and prone to failure.


\secref{theory} develops the theoretical framework for the globalization of
successive linearization methods with a differentiable merit function and 
shows that the resulting algorithm is well-defined for arbitrary 
complementarity problems, has strong global convergence properties, and is 
locally fast convergent under a strong regularity assumption.  
\secref{implementation} then discusses the PATH 4.x implementation of this 
framework which uses the successive linearization method of 
\cite{ralph:global,robinson:newtons} to generate the directions.
\secref{path-results} compares the reliability of PATH 4.x to PATH 2.9
on several test sets, demonstrating that PATH 4.x has better numerical 
behavior.  Finally, \secref{user} documents the PATH 4.x code from
the perspective of a user.

\section{Algorithm and Theory}\label{sec:theory}
The first component of the globalization strategy developed is a reformulation
of the mixed complementarity problem using the Fischer-Burmeister function 
\cite{fischer:special}.  Recall from \secref{refm} that the Fischer-Burmeister
function $ \phi_{FB}(a,b) := \sqrt{a^2 + b^2} - a - b$ is an NCP-function.
To reformulate the mixed complementarity problem using this function we 
introduced a partitioning of the index set $I = \{1, \ldots, n\}$:
\begin{eqnarray*}
   I_l & := & \{ i \in I \, | \, - \infty < l_i < u_i = + \infty \}, \\
   I_u & := & \{ i \in I \, | \, - \infty = l_i < u_i < + \infty \}, \\
   I_{lu} & := & \{ i \in I \, | \, - \infty < l_i < u_i < + \infty \}, \\
   I_f & := & \{ i \in I \, | \, - \infty = l_i < u_i = + \infty \},
\end{eqnarray*}
where $ I_l, \, I_u, \, I_{lu} $ and $ I_f $ denote the set of 
indices $ i \in I $ where there are finite lower bounds only, finite
upper bounds only, finite lower and upper bounds and no finite bounds
on the variable $ x_i $, respectively.  Following the idea in 
\cite{billups:algorithms}, we then defined the operator 
$ \Phi : \rn \to \rn $ componentwise as:
\[
   \Phi_i (x) := \left\{ \begin{array}{ll}
   \phi_{FB} (x_i- l_i, F_i (x)) & \mbox{if } i \in I_l, \\
   - \phi_{FB} (u_i-x_i, - F_i (x)) & \mbox{if } i \in I_u, \\
   \phi_{FB} (x_i - l_i, \phi_{FB} (u_i - x_i, - F_i (x))) &
   \mbox{if } i \in I_{lu}, \\
   - F_i(x) & \mbox{if } i \in I_f.
   \end{array} \right.
\]
We then have the following characterization of solutions to the mixed
complementarity problem:
\begin{proposition}[\cite{billups:algorithms}]\label{MCPEqual}
Let $x^* \in \rn$ be given.  Then the following are equivalent:
\begin{itemize}
\item[(a)] $ x^* $ is a solution of the mixed complementarity problem.
\item[(b)] $\Phi(x^*) = 0 $.
\end{itemize}
\end{proposition}

\noindent
\secref{reform} contains information related to $\Phi$ and its
subdifferential and in particular shows that a strongly regular 
solution to the complementarity problem implies that $\Phi$ is 
BD-regular at that point.  This fact will be needed to establish
local convergence of the method proposed.

The corresponding merit function for this nonsmooth system of equations:
\[
   \Psi (x) := \frac{1}{2} \Phi (x)\trt \Phi (x) = 
   \frac{1}{2} \| \Phi (x) \|^2
\]
is continuously differentiable.  In the special case of nonlinear 
complementarity problems, both $ \Phi $ and $ \Psi $ are used to design 
unconstrained algorithms for the solution of the problem
\cite{deluca.facchinei.ea:semismooth,facchinei.soares:merit}.  Unfortunately, 
in many practical situations, the imposed bounds, $l$ and $u$, on the 
variables are important not only for the problem definition but also because 
the complementarity function $ F $ may not be defined outside of these bounds.
For example, applications that include fractional powers can cause severe 
difficulties if the function is evaluated outside the feasible region.

Therefore, the basic algorithmic framework developed does not consider 
$\Phi$ directly, but instead attempts to solve the bound constrained 
optimization reformulation
$$
   \min \, \Psi (x) \quad \mbox{s.t.} \quad x \in [\ell,u].
$$
\secref{smooth} shows that a constrained stationary point of this
problem is a solution of the mixed complementarity problem under exactly 
the same assumptions used to prove a similar result for unconstrained 
stationary points of $\Psi$.

\secref{framework} then presents the algorithmic framework which only assumes 
we have a pre-existing feasible method that is locally well-defined and 
superlinearly convergent.  The algorithmic framework generates iterates that 
lie within the bounds, resorting to a projected gradient step for the bound 
constrained problem whenever the pre-existing method fails to provide 
a direction satisfying a sufficient decrease condition.  We prove that this 
method is well-defined as well as globally and locally fast convergent.

\subsection{Equation Properties}\label{sec:reform}
The function $ \Phi $ is not differentiable everywhere. However, it
is locally Lipschitzian and therefore has a nonempty generalized
Jacobian in the sense of Clarke \cite{clarke:optimization}.
We will use the following overestimation of this generalized Jacobian:
\begin{proposition}[\cite{billups:algorithms}]\label{GenJac}
We have
$$
   \partial \Phi (x)\trt \subseteq \{ D_a (x) + \nabla F(x) D_b (x) \},
$$
where $ D_a (x) \in \rntn $ and $ D_b (x) \in \rntn $ are diagonal
matrices whose diagonal elements are defined as follows:
\begin{itemize}
   \item[(a)] If $ i \in I_l $, then if $ (x_i-l_i, F_i (x)) \not = (0,0) $, 
      \begin{eqnarray*}
         (D_a)_{ii} (x) & = & \frac{x_i-l_i}{\| (x_i-l_i, F_i (x)) \|}-1, \\
         (D_b)_{ii} (x) & = & \frac{F_i(x)}{\| (x_i-l_i, F_i (x)) \|}-1
      \end{eqnarray*}
      but if $ (x_i-l_i, F_i (x)) = (0,0) $,
      $$
         ((D_a)_{ii} (x), (D_b)_{ii} (x)) \in \{ (\xi -1, \rho -1) 
         \in \Re^2 | \, \| (\xi, \rho) \| \le 1 \} .
      $$
   \item[(b)] If $ i \in I_u $, then
      if $ (u_i-x_i, -F_i (x)) \not = (0,0) $,
      \begin{eqnarray*}
         (D_a)_{ii} (x) & = & \frac{u_i-x_i}{\| (u_i-x_i, -F_i (x)) \|}-1, \\
         (D_b)_{ii} (x) & = & \frac{-F_i(x)}{\| (u_i-x_i, -F_i (x)) \|}-1
      \end{eqnarray*}
      but if $ (u_i-x_i, -F_i (x)) = (0,0) $,
      $$
         ((D_a)_{ii} (x), (D_b)_{ii} (x)) \in \{ (\xi -1, \rho -1) 
         \in \Re^2 | \, \| (\xi, \rho) \| \le 1 \} .
      $$ 
   \item[(c)] If $ i \in I_{lu} $, then
      $$
         (D_a)_{ii}(x) = a_i (x)+b_i (x) c_i (x), \qquad
         (D_b)_{ii}(x) = b_i (x) d_i (x) .
      $$
      Here,
      if $ (x_i-l_i, \phi (u_i-x_i, -F_i(x))) \not= (0,0) $,
      \begin{eqnarray*}
         a_i (x) & = & \frac{x_i-l_i}{\| (x_i-l_i, \phi (u_i-x_i, -F_i(x))) 
         \|} -1, \\
         b_i (x) & = & \frac{\phi (u_i-x_i, -F_i(x))}{\| 
         (x_i-l_i, \phi (u_i-x_i, -F_i(x))) \|} -1
      \end{eqnarray*}
      but if $ (x_i-l_i, \phi (u_i-x_i, -F_i(x))) = (0,0) $, 
      $$
         (a_i(x), b_i(x)) \in \{ (\xi -1, \rho -1) \in \Re^2 | \,
         \| (\xi, \rho) \| \le 1 \} .
      $$
      Further, if $ (u_i-x_i, -F_i (x)) \not= (0,0) $, then
      $$
         c_i (x) = \frac{x_i - u_i}{\| (u_i-x_i, -F_i (x)) \|}+1, \quad
         d_i (x) = \frac{F_i (x)}{\| (u_i-x_i, -F_i (x)) \|}+1
      $$
      but if $ (u_i-x_i, -F_i (x)) = (0,0) $,
      $$
         (c_i(x), d_i(x)) \in \{ (\xi +1, \rho +1) \in \Re^2 | \,
         \| (\xi, \rho) \| \le 1 \} .
      $$
   \item[(d)] If $ i \in I_f $, then
      $
         (D_a)_{ii} (x) = 0, \, (D_b)_{ii} (x) = -1
      $. 
\end{itemize}
\end{proposition}

\noindent
The statement of Proposition \ref{GenJac} is rather lengthy
because we have to take into account the definition of $ \Phi $
using the four different index sets $ I_l, I_u, I_{lu} $ and $ I_f $.
However, Proposition \ref{GenJac} is extremely important in the
subsequent analysis and will be used several times within the proofs of
some important results established in this and the next section. 

The remainder of this section is devoted to proving Theorem
\ref{BDregularity} which says that if $x^*$ is a strongly regular
solution to the mixed complementarity problem, then $\Phi$ is
BD-regular at $x^*$.  The motivation for this result is to 
establish the local convergence of the algorithm given 
in \secref{framework}.  To that end, let $ x^* \in \rn $ be a solution of the 
mixed complementarity problem and partition $I$ as follows:
\begin{eqnarray*}
   \alpha(x^*) & := & \{ i \, | \, l_i < x_i^* < u_i, F_i (x^*) = 0 \}, \\
   \beta(x^*) & := & \{ i \, | \, x_i^* \in \{ l_i, u_i \}, F_i (x^*) = 0 \}, \\
   \gamma(X^*) & := & \{ i \, | \, x_i^* \in \{ l_i, u_i \}, F_i (x^*) \not= 0 \}.
\end{eqnarray*}
The following result is then a simple consequence of Proposition \ref{GenJac}. 

\begin{lemma}\label{PropsH}
Let $ x^* \in \rn $ be a solution of the mixed 
complementarity problem. Furthermore, let $ H \in \partial \Phi (x^*) $
be any fixed matrix, $ H=D_a(x^*)+D_b(x^*) F'(x^*) $ with diagonal
matrices $ D_a (x^*) $ and $ D_b (x^*) $ as specified in Proposition
\ref{GenJac}. Then these diagonal matrices have the following properties:
\begin{itemize}
   \item[(a)] $ (D_a)_{ii} (x^*) = 0 $ and
      $ (D_b)_{ii} (x^*) = -1 $ for all $ i \in \alpha(x^*) $.
   \item[(b)] $ (D_a)_{ii} (x^*) \le 0, (D_b)_{ii} (x^*) \le 0, $ and
      $ (D_a)_{ii} (x^*) + (D_b)_{ii} (x^*) < 0 $ for all $ i \in \beta(x^*) $.
   \item[(c)] $ (D_a)_{ii} (x^*) = -1 $ and
      $ (D_b)_{ii} (x^*) = 0 $ for all $ i \in \gamma(x^*) $.
\end{itemize}
\end{lemma}

\noindent
{\bf Proof.}
If $ i \in \alpha(x^*) $, then we immediately obtain statement (a)
from Proposition \ref{GenJac} by considering the four possible cases
$ i \in I_l, i \in I_u, i \in I_{lu} $ and $ i \in I_f $ separately.

Next consider statement (c), i.e., assume that
$ i \in \gamma(x^*) $. Then we either have $ x_i^* = l_i $ and
$ F_i(x^*) > 0 $ or we have $ x_i^* = u_i $ and $ F_i (x^*) < 0 $.

First assume that $ x_i^*=l_i $ and $ F_i (x^*) > 0 $. 
Then the index $ i $ necessarily belongs to $ I_l $ or to $ I_{lu} $.
If $ i \in I_l $, we obtain from Proposition \ref{GenJac} that
$ (D_a)_{ii} (x^*)=-1 $ and $ (D_b)_{ii} (x^*) = 0 $. On the other hand,
if $ i \in I_{lu} $, we get from Proposition \ref{GenJac}, together with
the observation that $ \phi (a,b)>0 $ outside the nonnegative orthant, that
$$
   (D_a)_{ii} (x^*) = a_i (x^*)+b_i(x^*) c_i(x^*) = -1 + 0 \cdot c_i (x^*) = -1
$$
and
$$
   (D_b)_{ii} (x^*) = b_i (x^*) d_i (x^*) = 0 \cdot d_i (x^*) = 0.
$$

The case where $ x_i^*= u_i $ and $ F_i (x^*)< 0 $ can be proven in 
a similar manner.  Furthermore, statement (b) also follows using an 
identical argument.
\hfill$\Box$ \\

\noindent
We next restate a useful characterization of the strong regularity condition 
\cite{robinson:strongly} in the context of mixed complementarity problems.
The strong regularity condition is essentially a generalization of the
nonsingularity of the Jacobian assumption used for nonlinear equations.
In the case of a nonlinear complementarity problem, 
MCP($F,\{0\}^n,\{\infty\}^n$), this characterization reduces to the 
standard characterization in \cite{robinson:strongly}.

\begin{proposition}[\cite{facchinei.fischer.ea:semismooth}]\label{Charregularity}
The following are equivalent:
\begin{itemize}
   \item[(a)] $ x^* $ is a strongly regular solution of the mixed
      complementarity problem.
   \item[(b)] The submatrix $ F'(x^*)_{\alpha \alpha} $ is nonsingular,
      and the Schur-complement
      $$
         F'(x^*)_{\alpha \cup \beta, \alpha \cup \beta} (x^*)/
         F'(x^*)_{\alpha \alpha} := F'(x^*)_{\beta \beta} -
         F'(x^*)_{\beta \alpha} F'(x^*)_{\alpha \alpha}^{-1}
         F'(x^*)_{\alpha \beta}
      $$
      is a $ P $-matrix.
\end{itemize}
\end{proposition}

\noindent
In order to establish a nonsingularity result for the generalized Jacobian
$ \partial \Phi (x^*) $ at a strongly regular solution of the mixed
complementarity problem, we also need the following result; see 
\cite{gowda.sznajder:generalizations} for several extensions.

\begin{proposition}[\cite{kanzow.kleinmichel:new}]\label{CharP}
Let $M \in \real^{n \times n}$ be given.  The following are equivalent:
\begin{itemize}
\item[(a)] $M$ is a P-matrix.
\item[(b)] For all negative semidefinite diagonal matrices 
$D_a, D_b \in \Re^{n \times n} $ with $ D_a+D_b $ is negative definite:
\[
   D_a + D_b M
\]
is nonsingular.
\end{itemize}
\end{proposition}

\noindent
We are now able to prove the main result of this section.

\begin{theorem}\label{BDregularity}
If $ x^* $ is a strongly regular solution of the mixed complementarity problem,
then all elements $ H \in \partial \Phi (x^*) $ are nonsingular.  In
particular, $\Phi$ is BD-regular at $x^*$.
\end{theorem}

\noindent
{\bf Proof.}
Let $ H \in \partial \Phi (x^*) $. Then, by Proposition \ref{GenJac}, there
exist diagonal matrices $ D_a (x^*) , D_b (x^*) \in \rntn $ 
such that
\begin{equation}\label{Repres}
   H= D_a (x^*)+D_b(x^*) F'(x^*).
\end{equation}
Hence, if we write
$$
   D_a (x^*) = \left(\begin{array}{ccc}
   (D_a)_{\alpha \alpha} (x^*) & 0 & 0 \\
   0 & (D_a)_{\beta \beta} (x^*) & 0 \\
   0 & 0 & (D_a)_{\gamma \gamma} (x^*)
   \end{array} \right), 
$$
$$
   D_b (x^*) = \left(\begin{array}{ccc}
   (D_b)_{\alpha \alpha} (x^*) & 0 & 0 \\
   0 & (D_b)_{\beta \beta} (x^*) & 0 \\
   0 & 0 & (D_b)_{\gamma \gamma} (x^*)
   \end{array} \right) 
$$
and
$$
   F'(x^*) = \left(\begin{array}{ccc}
   F'(x^*)_{\alpha \alpha} &  F'(x^*)_{\alpha \beta} &
   F'(x^*)_{\alpha \gamma} \\
   F'(x^*)_{\beta \alpha} & F'(x^*)_{\beta \beta} &
   F'(x^*)_{\beta \gamma} \\
   F'(x^*)_{\gamma \alpha} & F'(x^*)_{\gamma \beta} &
   F'(x^*)_{\gamma \gamma}
   \end{array} \right) 
$$
and if we take into account Lemma \ref{PropsH}, the homogeneous linear
system $ Hd=0 $ can be rewritten as
\begin{eqnarray}
   F'(x^*)_{\alpha \alpha} d_{\alpha} + 
   F'(x^*)_{\alpha \beta} d_{\beta} +
   F'(x^*)_{\alpha \gamma} d_{\gamma} & = & 0_{\alpha}, \label{HS1} \\
   (D_a)_{\beta \beta} (x^*) d_{\beta} + (D_b)_{\beta \beta} (x^*) \left[
   F'(x^*)_{\beta \alpha} d_{\alpha} + 
   F'(x^*)_{\beta \beta} d_{\beta} +
   F'(x^*)_{\beta \gamma} d_{\gamma} \right] & = & 0_{\beta}, \label{HS2} \\
   - d_{\gamma} & = & 0_{\gamma}. \label{HS3}
\end{eqnarray}
Since $ d_{\gamma}=0 $ by (\ref{HS3}) and 
$ F'(x^*)_{\alpha \alpha} $ is nonsingular by assumption and
Proposition \ref{Charregularity}, we obtain from (\ref{HS1}):
\begin{equation}\label{HS4}
   d_{\alpha} = - F'(x^*)_{\alpha \alpha}^{-1} F'(x^*)_{\alpha \beta}
   d_{\beta}.
\end{equation}
Substituting (\ref{HS3}) and (\ref{HS4}) into (\ref{HS2}), we obtain
after some rearrangements:
\begin{equation}\label{HS5}
   \left[ (D_a)_{\beta \beta} (x^*)+ (D_b)_{\beta \beta} (x^*)
   (F'(x^*)_{\alpha \cup \beta, \alpha \cup \beta}/F'(x^*)_{\alpha \alpha})
   \right] d_{\beta} = 0_{\beta}.
\end{equation}
Since the Schur complement 
$ F'(x^*)_{\alpha \cup \beta, \alpha \cup \beta}/F'(x^*)_{\alpha \alpha} $
is a $ P $-matrix by assumption and Proposition \ref{Charregularity}
and since, by Lemma \ref{PropsH} (b), 
the diagonal matrices $ (D_a)_{\beta \beta} (x^*) $
and $ (D_b)_{\beta \beta} (x^*) $ are negative semidefinite with a
negative definite sum, it follows from Proposition \ref{CharP} that the
coefficient matrix in (\ref{HS5}) is nonsingular. Hence we obtain
$ d_{\beta} = 0_{\beta} $. This, in turn, implies $ d_{\alpha} = 0_{\alpha} $
by (\ref{HS4}). Reference to (\ref{HS3}) shows that $ d = 0 $, so that
$ H $ is nonsingular.  Since 
$\partial_B \Phi(x^*) \subseteq \partial \Phi(x^*)$ it follows that 
$\Phi$ is BD-regular at $x^*$.
\hfill$\Box$ \\

\subsection{Merit Function Properties}\label{sec:smooth}
We now investigate properties of the residual merit function
$$
   \Psi (x) = \frac{1}{2} \Phi (x)\trt \Phi (x)
$$
associated with the equation operator $ \Phi $.  Despite the fact that 
$ \Phi $ is nondifferentiable in general, the merit function $ \Psi $ is 
continuously differentiable everywhere. 
\begin{proposition}[\cite{billups:algorithms}]\label{PsiDiff}
The function $ \Psi $ is continuously differentiable with gradient
$ \nabla \Psi (x) = H\trt \Phi (x) $ for an arbitrary
$ H \in \partial \Phi (x) $.
\end{proposition}

\noindent
We next provide a stationary point result for the unconstrained
reformulation
$$
   \min \, \Psi (x), \qquad x \in \rn,
$$
of the mixed complementarity problem. 
To this end, we need the following characterization of $ P_0 $-matrices.
Generalizations of this result can be found in
\cite{chen.harker:non-interior,gowda.sznajder:generalizations}.

\begin{proposition}\label{CharP0}
Let $M \in \real^{n \times n}$ be given.  The following are equivalent:
\begin{itemize}
\item[(a)] $M$ is a $P_0$-matrix.
\item[(b)] For all negative definite diagonal matrices 
$D_a, D_b \in \Re^{n \times n} $:
\[
   D_a + D_b M
\]
is nonsingular.
\end{itemize}
\end{proposition}

\noindent
We will also need the following technical result.

\begin{lemma}\label{Tech2}
Let $ x \in \rn $ be arbitrary and $ H \in \partial \Phi (x), 
H= D_a (x) + D_b(x) F'(x) $ with diagonal matrices
$ D_a (x), D_b (x) \in \rntn $ as defined in Proposition \ref{GenJac}. Then 
the following two statements hold:
\begin{itemize}
   \item[(a)] For all $ i \in I $,
        $ [D_a (x) \Phi (x)]_i [D_b (x) \Phi (x)]_i \ge 0 $ .
   \item[(b)] For each $ i \not\in I_f $,
        $ [D_a (x) \Phi (x)]_i =0 \Longleftrightarrow
      [D_b (x) \Phi (x)]_i =0 \Longleftrightarrow \Phi_i (x)=0 $ .
\end{itemize}
\end{lemma}

\noindent
{\bf Proof.}
(a) By considering the four possible cases $ i \in I_l, i \in I_u, i \in
I_{lu} $ and $ i \in I_f $, it is easy to see that $ (D_a)_{ii} (x) \le 0 $
and $ (D_b)_{ii} (x) \le 0 $ for all $ i \in I $. Hence
$$
   [D_a (x) \Phi (x)]_i [D_b (x) \Phi (x)]_i =
   (D_a)_{ii} (x) (D_b)_{ii} (x) \Phi_i (x)^2 \ge 0
$$
for all $ i \in I $. \\ \smallskip

\noindent
(b) If $ \Phi_i (x) = 0 $, we immediately have
$$
   [D_a (x) \Phi (x)]_i =0 \quad \mbox{and} \quad
   [D_b (x) \Phi (x)]_i =0.
$$
Conversely, assume that $ [D_a (x) \Phi (x)]_i =0 $ for some index
$ i \not\in I_f $ (the proof is analogous if $ [D_b (x) \Phi (x)]_i =0 $).
Then
$$
    (D_a)_{ii} (x)=0 \quad \mbox{or} \quad \Phi_i (x) =0.
$$
In the latter case, there is nothing to show. So suppose that
$ (D_a)_{ii} (x) = 0 $. Due to the definition of $ D_a (x) $, we distinguish
three cases. \\ \smallskip

\noindent
{\em Case 1:} $ i \in I_l $. \\
If $ (x_i-l_i, F_i(x))=(0,0) $, then $ \Phi_i (x)=0 $ follows immediately
from the definition of the operator $ \Phi $. Otherwise Proposition
\ref{GenJac} (a) gives
$$
   0 = (D_a)_{ii} (x) = \frac{x_i-l_i}{\| (x_i-l_i, F_i(x)) \|}-1.
$$
This implies $ x_i-l_i > 0 $ and $ F_i(x)=0 $, so that $ \Phi_i(x)=0 $
in view of the very definition of $ \Phi $ and the NCP-property of
the function $ \phi $. \\ \smallskip

\noindent
{\em Case 2:} $ i \in I_u $. \\
The proof of this case is very similar to the one given for Case 1 and
we therefore omit the details. \\ \smallskip

\noindent
{\em Case 3:} $ i \in I_{lu} $. \\
If $ (x_i-l_i, \phi (u_i-x_i, -F_i(x)) = (0,0) $, we are done.
So assume that $ (x_i-l_i, \phi (u_i-x_i, -F_i(x)) \not= (0,0) $.
Then Proposition \ref{GenJac} (c) gives
\begin{equation}\label{Daii}
   0 = (D_a)_{ii} (x) = a_i(x)+b_i(x) c_i(x)
\end{equation}
with certain numbers $ a_i(x), b_i(x) $ and $ c_i(x) $ specified in
Proposition \ref{GenJac} (c). Since it follows immediately from this
Proposition that
$$
   a_i(x) \le 0, \enspace b_i(x) \le 0 \quad \mbox{and} \quad
   c_i(x) \ge 0,
$$
the right-hand side of (\ref{Daii}) is the sum of two nonpositive 
expressions which can therefore be equal to zero only if
$
   a_i(x) = 0
$.
This, however, implies
$x_i-l_i>0$ and $\phi (u_i-x_i, -F_i(x)) = 0$
in view of the definition of $ a_i(x) $ given in Proposition \ref{GenJac} (c).
Hence $ \Phi_i(x)=0 $ by the NCP-property of the function $ \phi $.
\hfill$\Box$\\

\noindent
Note that Lemma \ref{Tech2} (b) does not hold for indices $ i \in I_f $
since $ (D_b)_{ii}(x)=-1 $ for all $ i \in I_f $ in view of
Proposition \ref{GenJac} (d).

Proposition \ref{CharP0} and Lemma \ref{Tech2} are used to prove the first 
major result of this section. In this result and in the remainder of this 
section, we use the short-hand notation $ \nabla F(x^*)_{ff} $ to denote the 
submatrix $ \nabla F(x^*)_{I_f I_f} $. A similar notation is used for 
submatrices and subvectors defined by other index sets.

\begin{theorem}\label{UStat}
Let $ x^* \in \rn $ be a stationary point of $ \Psi $. Assume that
\begin{itemize}
   \item[(a)] the principal submatrix $ \nabla F(x^*)_{ff} $ is nonsingular,
      and
   \item[(b)] the Schur complement $ \nabla F(x^*)/\nabla F(x^*)_{ff} $
      is a $ P_0 $-matrix.
\end{itemize}
Then $ x^* $ is a solution of the mixed complementarity problem.
\end{theorem}

\noindent
{\bf Proof.}
Let $ x^* $ be a stationary point of $ \Psi $. Then, by Proposition
\ref{PsiDiff}, we have
\begin{equation}\label{S1}
   H\trt \Phi (x^*) = \nabla \Psi (x^*) = 0
\end{equation} 
for an arbitrary $ H \in \partial \Phi (x^*) $. By Proposition \ref{GenJac},
there exist diagonal matrices $ D_a (x^*), D_b (x^*) \in \rntn $ such that
$$
   H = D_a (x^*) + D_b (x^*) F'(x^*).
$$
Therefore, (\ref{S1}) can be written as
\begin{equation}\label{S2}
   \left[ D_a(x^*) + \nabla F(x^*) D_b (x^*) \right] \Phi (x^*) = 0.
\end{equation}
Writing
$$
   D_a (x^*) = \left(\begin{array}{cc}
   (D_a)_{ff}(x^*) & 0 \\
   0 & (D_a)_{\bar f \bar f} (x^*) 
   \end{array} \right),
$$
$$
   D_b (x^*) = \left(\begin{array}{cc}
   (D_b)_{ff}(x^*) & 0 \\
   0 & (D_b)_{\bar f \bar f} (x^*)
   \end{array} \right)
$$
and
$$
   \nabla F(x^*) = \left(\begin{array}{cc}
   \nabla F(x^*)_{ff} & \nabla F(x^*)_{f \bar f} \\
   \nabla F(x^*)_{\bar f f} & \nabla F(x^*)_{\bar f \bar f}
   \end{array} \right),
$$
where $ I_{\bar f} := I \setminus I_f $, and taking into account that
\begin{eqnarray*}
   (D_a)_{ii} (x^*) & = & 0 \qquad \forall i \in I_f, \\
   (D_b)_{ii} (x^*) & = & -1 \qquad \forall i \in I_f
\end{eqnarray*}
by Proposition \ref{GenJac}, we can rewrite (\ref{S2}) as
\begin{eqnarray}
   - \nabla F(x^*)_{ff} \Phi (x^*)_f + \nabla F(x^*)_{f \bar f}
   (D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f} & = & 0_f, \label{S3} \\
   (D_a)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f} - \nabla F(x^*)_{\bar f f}
   \Phi (x^*)_f + \nabla F(x^*)_{\bar f \bar f} (D_b)_{\bar f \bar f} (x^*)
   \Phi (x^*)_{\bar f} & = & 0_{\bar f}. \label{S4} 
\end{eqnarray}
Due to the assumed nonsingularity of $ \nabla F(x^*)_{ff} $, we obtain from
(\ref{S3}):
\begin{equation}\label{S5}
   \Phi (x^*)_f = \nabla F(x^*)_{ff}^{-1} \nabla F(x^*)_{f \bar f}
   (D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}.
\end{equation}
Substituting this expression into (\ref{S4}) and rearranging terms gives
\begin{equation}\label{S6}
   \left[ (D_a)_{\bar f \bar f} (x^*)+ (\nabla F(x^*)/\nabla F(x^*)_{ff})
   (D_b)_{\bar f \bar f} (x^*) \right] \Phi (x^*)_{\bar f} = 0_{\bar f}.
\end{equation}
Define the index sets:
\begin{eqnarray*}
\begin{array}{lcl}
f_0 & := & \{i \mid [(D_a)_{\bar f \bar f}(x^*)]_{ii} = 0 \text{ or } [(D_b)_{\bar f \bar f}(x^*)]_{ii} = 0\}\\
\bar{f_0} & := & \{i \mid [(D_a)_{\bar f \bar f}(x^*)]_{ii} < 0 \text{ and } [(D_b)_{\bar f \bar f}(x^*)]_{ii} < 0\}.
\end{array}
\end{eqnarray*}
Using Proposition \ref{GenJac}, we have that 
$f_0 \cup \bar{f_0} = \{1, \ldots, k\}$ where $k$ is the cardinality of 
$\bar f$.  Then, by Lemma \ref{Tech2}, $[\Phi(x^*)_{\bar f}]_{f_0} = 0$.  
Consequently, (\ref{S6}) reduces to the system:
\begin{equation}\label{S6a}
   \left[ [(D_a)_{\bar f \bar f}(x^*)]_{\bar{f_0} \bar{f_0}} + [\nabla F(x^*)/\nabla F(x^*)_{ff}]_{\bar{f_0} \bar{f_0}}
   [(D_b)_{\bar f \bar f}(x^*)]_{\bar{f_0} \bar{f_0}} \right] 
\Phi (x^*)_{\bar{f_0}} = 0_{\bar{f_0}}.
\end{equation}
Since $(\nabla F(x^*)/\nabla F(x^*)_{ff})$ is a $P_0$-matrix by 
assumption, it follows that the principle submatrix, 
$[\nabla F(x^*)/\nabla F(x^*)_{ff}]_{\bar{f_0} \bar{f_0}}$, is also
$P_0$-matrix.  Then by Proposition 
\ref{CharP0} the coefficient matrix in (\ref{S6a}) is nonsingular.  Hence, 
$[\Phi (x^*)_{\bar f}]_{\bar{f_0}} = 0$.  Therefore,
$ \Phi (x^*)_{\bar f} = 0_{\bar f} $. But then (\ref{S5}) implies 
$ \Phi (x^*)_f = 0_f $. Hence, $ \Phi (x^*) = 0 $, i.e., $ x^* $ is a 
solution of the mixed complementarity problem by Proposition \ref{MCPEqual}.
\hfill$\Box$ \\

\noindent
We now provide a sufficient condition for a stationary point of
the constrained reformulation of the mixed complementarity problem
\begin{equation}\label{CR}
   \min \, \Psi (x) \quad \mbox{s.t.} \quad x \in [\ell,u]
\end{equation}
to be a global minimum.  In fact, this result is of much more importance than
the unconstrained stationary point result given in Theorem \ref{UStat}.
However, Theorem \ref{UStat} will be used to establish the constrained 
stationary point result.

\begin{theorem}\label{CStat}
Let $ x^* \in \rn $ be a stationary point of the constrained
reformulation (\ref{CR}) of the mixed complementarity problem. Assume that
\begin{itemize}
   \item[(a)] the principal submatrix $ \nabla F(x^*)_{ff} $ is nonsingular,
      and
   \item[(b)] the Schur complement $ \nabla F(x^*)/\nabla F(x^*)_{ff} $
      is a $ P_0 $-matrix.
\end{itemize}
Then $ x^* $ is a solution of the mixed complementarity problem.
\end{theorem}

\noindent
{\bf Proof.}
Since $ x^* $ is a stationary point of the reformulation (\ref{CR}),
it satisfies the following conditions (which themselves form a 
mixed complementarity problem):
\begin{equation}\label{rel}
\begin{array}{rcl}
   x_i^* = l_i & \Longrightarrow & [\nabla \Psi (x^*)]_i \ge 0, \\
   x_i^* = u_i & \Longrightarrow & [\nabla \Psi (x^*)]_i \le 0, \\
   x_i^* \in (l_i,u_i) & \Longrightarrow & [\nabla \Psi (x^*)]_i = 0.
\end{array}
\end{equation}
The main part in proving that $ x^* $ is already a solution of the
mixed complementarity problem consists in showing that we actually
have $ [\nabla \Psi (x^*)]_i = 0 $ for all $ i \in I $.

The proof is by contradiction, so assume that $ \nabla \Psi (x^*) \neq 0 $.
Since $ \nabla \Psi (x^*) $ can be written as
$$
   \nabla \Psi (x^*) = H\trt \Phi (x^*) =
   D_a(x^*) \Phi (x^*) + \nabla F(x^*) D_b(x^*) \Phi (x^*)
$$
for a matrix $ H \in \partial \Phi (x^*) $ and certain diagonal matrices
$ D_a (x^*), D_b (x^*) \in \rntn $ by Propositions \ref{GenJac} and
\ref{PsiDiff}, and since we necessarily have 
$$
   [ \nabla \Psi (x^*) ]_f = 0_f
$$
because of (\ref{rel}), we can follow the argument used in the proof
of Theorem \ref{UStat} in order to show that
$$
   \Phi (x^*)_f = \nabla F(x^*)_{ff}^{-1} \nabla F(x^*)_{f \bar f}
   (D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f} .
$$
Substituting this into the expression for $ [\nabla \Psi (x^*)
]_{\bar f} $ and rearranging terms, we obtain
\begin{equation}\label{express}
   [\nabla \Psi (x^*)]_{\bar f} = (D_a)_{\bar f \bar f} (x^*) 
   \Phi (x^*)_{\bar f}
   + (\nabla F(x^*)/ \nabla F(x^*)_{ff}) (D_b)_{\bar f \bar f} (x^*)
   \Phi (x^*)_{\bar f} .
\end{equation}
Since $ \nabla \Psi (x^*) \neq 0 $ by assumption, there exists an index 
$ i \in I_{\bar f} $ such that either
\begin{equation}\label{E1}
   x_i^* = l_i \quad \mbox{and} \quad [\nabla \Psi (x^*)]_i > 0
\end{equation}
or
\begin{equation}\label{E2}
   x_i^* = u_i \quad \mbox{and} \quad [\nabla \Psi (x^*)]_i < 0.
\end{equation}
Now it follows easily from Proposition \ref{GenJac} that if $
x_i^*=l_i $, then
$$
   [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_i \le 0
$$
and that if $ x_i^* =u_i $,
$$
   [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_i \ge 0.
$$
Therefore, if we premultiply $ [\nabla \Psi (x^*)]_i $ in
(\ref{E1}) and (\ref{E2}) by 
$ [(D_b)_{\bar f \bar f} \Phi (x^*)_{\bar f}]_i $ 
and substitute the $i$th component from the
expression (\ref{express}) 
for $ [\nabla \Psi (x^*)]_i $, we obtain in both cases that
\begin{multline}\label{E5}
   [(D_a)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_i 
   [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_i + \\
   [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_i 
   [(\nabla F(x^*)/\nabla F(x^*)_{ff}) (D_b)_{\bar f \bar f} (x^*) 
   \Phi (x^*)_{\bar f}]_i \le 0.
\end{multline}
Note that this inequality holds for all indices $ i \in I_{\bar f} $ such that
$ [\nabla \Psi (x^*)]_i \neq 0 $. In addition, we can show in a similar way 
that the equality
\begin{multline}\label{E6}
   [(D_a)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_i 
   [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_i + \\
   [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_i 
   [(\nabla F(x^*)/\nabla F(x^*)_{ff}) (D_b)_{\bar f \bar f} (x^*) 
   \Phi (x^*)_{\bar f}]_i = 0
\end{multline}
holds for all indices $ i \in I_{\bar f} $ with $ [\nabla \Psi (x^*)]_i = 0 $.

Since $ \nabla \Psi (x^*) \neq 0 $ by assumption
and since we already know that $ [\nabla \Psi (x^*)]_f = 0 $, it follows
immediately from (\ref{express}) and Lemma \ref{Tech2} (b) that
$ (D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f} $ is a nonzero vector.
Therefore, since the Schur complement $ \nabla F(x^*)/\nabla F(x^*)_{ff} $
is a $ P_0 $-matrix by assumption, there
exists an index $ i_0 \in I_{\bar f} $ such that
\begin{multline}\label{contra}
   [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_{i_0} \neq 0 
   \quad \mbox{and} \quad \\
   [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_{i_0} 
   [(\nabla F(x^*)/\nabla F(x^*)_{ff}) (D_b)_{\bar f \bar f} (x^*) 
   \Phi (x^*)_{\bar f}]_{i_0} \ge 0.
\end{multline}
Now Lemma \ref{Tech2} (a), (\ref{E5}), (\ref{E6}) 
and (\ref{contra}) imply that
$$
   0 = [(D_a)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_{i_0} 
   [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_{i_0}
$$
and therefore $ [(D_b)_{\bar f \bar f} (x^*) \Phi (x^*)_{\bar f}]_{i_0} = 0 $ 
by Lemma \ref{Tech2} (b).
This, however, contradicts the choice of the index $ i_0 $ in (\ref{contra}).

Hence we must have
$ \nabla \Psi (x^*) =0 $, so that Theorem
\ref{UStat} gives the desired result that $ x^* $ is a solution of the
mixed complementarity problem.
\hfill$\Box$\\

\noindent
We note that, if we apply the main results of this section to the standard
nonlinear complementarity problem, then we obtain some known properties
\cite{deluca.facchinei.ea:semismooth,facchinei.soares:merit,fischer:constrained} of the merit function $ \Psi $.

\subsection{Algorithmic Framework}\label{sec:framework}
We now present the algorithm used to solve the mixed complementarity problem 
and corresponding global and local convergence theory.  We assume that we 
have a basic algorithm, which we will call {\bf Algorithm LM} (for local 
method), with the following two properties:
\begin{itemize}
   \item[(a)] Given any point $ x^k \in [\ell,u] $, if Algorithm LM is able to 
      compute a search direction $ d^k \in \rn $, then this direction 
      satisfies $ x^k + d^k \in [\ell,u] $;
   \item[(b)] Given any sequence $ \{ x^k \} $ converging to a strongly
      regular solution $ x^* $ of the mixed complementarity problem, 
      Algorithm LM is able to compute a search direction $ d^k $ for all
      $ x^k $ sufficiently close to $ x^* $, and this direction has the 
      property that $ \| x^k + d^k - x^* \| = o (\| x^k - x^* \|) $.
\end{itemize}
Property (a) is a very weak assumption; it does not even assume that
Algorithm LM is able to do anything at an arbitrary given point $ x^k $.
For example, the subproblem might be inconsistent.
However, if Algorithm LM is able to compute a search direction, we assume 
that it computes a search direction such that, if we take the full step, 
the new point $ x^k + d^k $ stays in the feasible set $ [\ell,u] $. Note that,
due to the convexity of the feasible set $ [\ell,u] $, this guarantees that 
all the points $ x^k + t_k d^k, t_k \in [0,1] $ are feasible, too. On the 
other hand, property (b) states that, under the strong regularity condition, 
Algorithm LM is locally well-defined and generates a locally superlinearly
convergent search direction. Hence we view Algorithm LM as a feasible and 
locally superlinearly convergent method for the solution of the mixed
complementarity problem.  Note, however, that the above two conditions say 
nothing about the way in which we generate the sequence $ \{ x^k \} $.

Several methods satisfy the above two conditions.  For example, one may take 
the Josephy-Newton method \cite{josephy:newtons,mathiesen:algorithm}
or the alternative method suggested in \cite{robinson:newtons,ralph:global}.
The NE/SQP method \cite{pang.gabriel:nesqp} is another possible candidate 
as is the inexact QP-based solver in \cite{kanzow:inexact}.  These latter two 
methods have been used to solve the standard complementarity problem only, 
but it is not difficult to extend both methods to mixed complementarity 
problems \cite{billups:algorithms,billups.ferris:qpcomp}.

The algorithm presented globalizes Algorithm LM using the merit function 
$ \Psi $ to measure any progress. If the point generated by Algorithm LM has 
a merit function value of $ \Psi $ sufficiently smaller than the previous 
one, it 
is accepted as the new iterate. Otherwise a projected gradient step for the 
smooth merit function $ \Psi $ is taken. In this way, we guarantee that all 
iterates stay in the feasible set $ [\ell,u] $. In effect, the class of 
methods is an algorithmic framework for the solution of the box constrained 
optimization problem
\begin{equation}\label{BCO}
   \min \, \Psi (x) \quad \mbox{s.t.} \quad x \in [\ell,u].
\end{equation}

We now provide a detailed statement of the algorithm, where the projection 
of an arbitrary point $ x \in \rn $ on the feasible set
$ [\ell,u] $ is denoted by $ \pi_{[\ell,u]}(x) $.

{\bf Algorithm GDM (General Descent Method)}
{\it
\begin{itemize}
   \item[(S.0)] (Initialization) \\
      Choose $ x^0 \in [\ell,u], s > 0, \beta \in (0,1), 
      \gamma \in (0,1) $ and set $ k:= 0 $.
   \item[(S.1)] (Termination Criterion) \\
      If $ x^k $ is a stationary point of (\ref{BCO}): STOP. 
   \item[(S.2)] (Compute Fast Search Direction) \\
      Use Algorithm LM to compute a search direction $ d^k $. If this is
      not possible or if the condition
      \begin{equation}\label{descent}
         \Psi (x^k+d^k) \le \gamma \Psi (x^k)
      \end{equation}
      is not satisfied, go to Step (S.4), else go to Step (S.3).
   \item[(S.3)] (Accept Fast Search Direction) \\
      Set $ x^{k+1}:=x^k+d^k, k \leftarrow k+1 $, and go to Step (S.1).
   \item[(S.4)] (Take Projected Gradient Step) \\
      Compute $ t_k = \max \{ s \beta^{i} | \, i = 0, 1, 2, \ldots \} $
      such that
      \begin{equation}\label{LSAlg}
         \Psi (x^k(t_k)) \le \Psi (x^k) - \sigma \nabla 
         \Psi (x^k)\trt (x^k - x^k (t_k)), 
      \end{equation}
      where $ x^k (t):= \pi_{[\ell,u]}(x^k - t \nabla \Psi (x^k)) $.
      Set $ x^{k+1}:= x^k (t_k), k \leftarrow k+1 $, and go to Step (S.1).
\end{itemize}
}

\noindent 
Other methods have attempted to use projected gradient steps in
conjunction with steps that give fast local convergence 
\cite{ferris.ralph:projected}.
Unfortunately, these hybrid algorithms are difficult to implement and
numerical testing has therefore only been carried out on small test examples.
A key difference in the approach outlined here is that the
implementation can be achieved as a modification to an existing
code. 

We now investigate the convergence properties of Algorithm GDM.
To this end, we always assume implicitly that Algorithm GDM does not 
terminate in a finite number of steps.  That is, none of the iterates 
$ x^k $ is a stationary point of (\ref{BCO}).

The global convergence analysis consists of two parts:
we first show that the algorithm is well-defined, and we then prove
that any accumulation point of a sequence $ \{ x^k \} $
generated by Algorithm GDM is a stationary point of the bound constrained 
optimization problem
\begin{equation}\label{COP}
   \min \, \Psi (x) \quad \mbox{s.t.} \quad x \in [\ell,u].
\end{equation}
Recall that Theorem \ref{CStat} gives a relatively mild condition for a 
stationary point of (\ref{COP}) to be a solution of the mixed complementarity 
problem.

\begin{theorem}\label{GlobConv}
Algorithm GDM is well-defined for an arbitrary mixed complementarity
problem with a continuously differentiable function $ F $ defined on
an open set containing the rectangle $ [\ell,u] $.  Furthermore, every 
accumulation point of a sequence $ \{ x^k \} $ generated by Algorithm GDM is 
a stationary point of (\ref{COP}). 
\end{theorem}

\noindent
{\bf Proof.}
To prove the algorithm is well-defined,
we only have to show that the projected gradient step can be carried out
at each iteration.  That is, we can always find a finite step length
$ t_k > 0 $ satisfying condition (\ref{LSAlg}). However, since we assume
that none of the iterates $ x^k $ is a stationary point of (\ref{BCO}),
this follows, e.g., from Proposition 2.3.3 (a) in Bertsekas 
\cite{bertsekas:nonlinear}.

For the second part of the theorem, let $ x^* $ be an accumulation point 
of the 
sequence $ \{ x^k \}, $ and assume that $ \{ x^k \}_K $ is a subsequence 
converging to $ x^* $.  Suppose there are infinitely many $ k \in K $ such 
that $ x^{k+1} $ is generated by using a projected gradient step for all 
these $ k $.  Since the iterates $ x^k $ belong to the feasible set 
$ [\ell,u] $ for all $ k \in \Ne $ and since the sequence $ \{ \Psi (x^k) \} $
is monotonically decreasing, it is not difficult to see that the 
proof of Proposition 2.3.3 (b) in Bertsekas \cite{bertsekas:nonlinear} can 
be adapted in a straightforward manner to establish that $ x^* $
is a stationary point of the constrained reformulation (\ref{COP}).

Hence we can assume without loss of generality that all iterates $ k \in K $
satisfy the descent condition (\ref{descent}). Due to the monotonic decrease 
of the sequence $ \{ \Psi (x^k) \} $, this implies that the entire sequence 
$ \{ \Psi (x^k) \} $ converges to $ 0 $. In particular, in view of the 
definition of the merit function, we see that the accumulation point $ x^* $ 
is a solution of the mixed complementarity problem and hence also a stationary
point of problem (\ref{COP}).
\hfill$\Box$ \\

We now show that Algorithm GDM is locally Q-superlinearly convergent under 
the strong regularity condition \cite{robinson:strongly}. The proof is in 
two parts: we first show that the entire sequence generated by Algorithm GDM 
converges to a solution $ x^* $ if this solution satisfies the strong 
regularity assumption, and then we determine the rate of convergence. 
The critical tool to establish that the sequence converges is the
following proposition.

\begin{proposition}[\cite{more.sorensen:computing}]\label{MSLemma}
Assume that $ x^* $ is an isolated accumulation point of a sequence
$ \{ x^k \} $ (not necessarily generated by Algorithm GDM) such that
$ \{ \| x^{k+1}-x^k \| \}_K \to 0 $ for any subsequence $ \{ x^k \}_K $
converging to $ x^* $. Then the whole sequence $ \{ x^k \} $ converges to
$ x^* $.
\end{proposition}

\noindent
The basic device used to show Q-superlinear convergence of Algorithm GDM
is to prove that eventually there are no projected gradient steps, so the 
method inherits the local convergence properties of the locally superlinearly 
convergent Algorithm LM used in Step (S.2).  To simplify the proof, we invoke 
the following proposition.
 
\begin{proposition}[\cite{facchinei.soares:merit,kanzow.qi:qp-free}]\label{KQProp}
Let $ G: \rn \to \rn $ be locally Lipschitzian, $ x^* \in \rn $ with
$ G(x^*)=0 $ be such that all elements in $ \partial G(x^*) $ are
nonsingular, and assume that there are two sequences $ \{ x^k \} \subseteq
\rn $ and $ \{ d^k \} \subseteq \rn $ (not necessarily generated by
Algorithm GDM) with $ \{ x^k \} \to x^* $ and
$ \| x^k + d^k -x^* \| = o(\| x^k - x^* \|) $. Then
$ \| G(x^k + d^k) \| = o (\| G(x^k) \|) $.
\end{proposition}

\noindent
We can now state and prove the main local convergence result for
Algorithm GDM.  The convergence rate established depends critically on 
the main result of \secref{reform}, namely Theorem \ref{BDregularity}. 

\begin{theorem}\label{ConvLocal}
Let $ \{ x^k \} \subseteq \rn $ be a sequence generated by Algorithm
GDM. Assume that this sequence has an accumulation point $ x^* $
which is a strongly regular solution of the mixed complementarity problem.
Then the entire sequence $ \{ x^k \} $ converges to this point, and the
rate of convergence is Q-superlinear. 
\end{theorem}

\noindent
{\bf Proof.}
To establish convergence,
we first note that a strongly regular solution is an isolated
solution of the mixed complementarity problem, see \cite{robinson:strongly}. 
Since Algorithm GDM generates a decreasing sequence $ \{ \Psi (x^k) \} $
and $ x^* $ is a solution of the mixed complementarity problem, the
entire sequence $ \{ \Psi (x^k) \} $ converges to zero. Hence every
accumulation point of the sequence $ \{ x^k \} $ must be a solution of
the mixed complementarity problem. Therefore, the assumed strong regularity of
$ x^* $ implies that $ x^* $ is an isolated accumulation point of the 
sequence $ \{ x^k \} $. 

Now let $ \{ x^k \}_K $ denote any subsequence converging to $ x^* $.
Assume first that we take a projected gradient step for all sufficiently
large $k \in K$.  Then we obtain, using the nonexpansive property of the 
projection operator:
\begin{equation}\label{JNStep}
\begin{array}{rcl}
   \| x^{k+1}-x^k \| & = & \| x^k (t_k) - x^k \| \\
   & = & \| [x^k - t_k \nabla \Psi (x^k)]_+ - x^k \| \\
   & = & \| [x^k - t_k \nabla \Psi (x^k)]_+ - [x^k]_+ \| \\
   & \le & \| t_k \nabla \Psi (x^k) \| \\
   & \le & s \| \nabla \Psi (x^k) \| \\
   & \to & 0
\end{array}
\end{equation}
since $ x^* $ solves the mixed complementarity problem
so that $ x^* $ is a global minimizer and hence a stationary
point of the merit function $ \Psi $.

On the other hand, if we calculate and accept the search direction $ d^k $ 
generated by Algorithm LM for infinitely many $ k \in K $, then the 
assumptions on Algorithm LM and the assumed strong regularity imply
that $ \{ d^k \} \to 0 $ on this infinite subsequence, so that the updating 
rules from Algorithm GDM show that
\begin{equation}\label{PGStep}
   \| x^{k+1}-x^k \| = \| d^k \| \to 0
\end{equation}
on this subsequence. 
Combining (\ref{PGStep}) and (\ref{JNStep}), we immediately obtain that
$$
   \{ \| x^{k+1}-x^k \| \}_K \to 0.
$$
Hence the assumptions of Proposition \ref{MSLemma} are satisfied, and
convergence follows from that result.

In order to establish the Q-superlinear rate, we note that
the strong regularity of the solution
$ x^* $ and Theorem \ref{BDregularity} show that all elements in
$ \partial \Phi (x^*) $ are nonsingular. Since, in view of the assumptions
about the search directions $ d^k $ generated by Algorithm LM, we have
$ \| x^k + d^k - x^* \| = o (\| x^k - x^* \|) $ for these search directions,
Proposition \ref{KQProp} implies that
$$
   \| \Phi (x^k + d^k) \| = o (\| \Phi (x^k) \|)
$$
and therefore
$$
   \Psi (x^k + d^k) = o (\Psi (x^k)).
$$
This shows that the descent condition
$$
   \Psi (x^k + d^k) = \gamma \Psi (x^k)
$$
is eventually
satisfied in Step (S.2) of Algorithm GDM, i.e., for all $ k \in \Ne $
sufficiently large, Algorithm GDM does not take any projected gradient
steps. Hence Algorithm GDM has the same local convergence properties
as Algorithm LM. Since $ x^{k+1}=x^k + d^k $, this means that $ \{ x^k \} $
converges Q-superlinearly to $ x^* $.
\hfill$\Box$\\

\noindent
Obviously, if the basic Algorithm LM is locally Q-quadratically convergent,
then the class of algorithms given in Algorithm GDM is also locally 
Q-quadratically convergent. Typically, this holds if we assume in addition 
that the Jacobian of $ F $ is locally Lipschitzian.

\section{Implementation Features}\label{sec:implementation}
PATH 4.x uses Algorithm GDM to globalize the Newton method from
\cite{ralph:global,robinson:newtons}.  This local method constructs
a point-based approximation of the normal map \cite{robinson:normal}:
\[ 
    F(\pi_{[\ell,u]}(x)) + x - \pi_{[\ell,u]}(x) 
\]
where $\pi_{[\ell,u]}(x)$ represents the projection of $x$ onto 
$[\ell,u]$.  The resulting linear complementarity problem is solved 
using a pivotal code related to Lemke's method \cite{lemke:complementary} 
to calculate a Newton point.

Two implementational points are of interest.  We first need to reconcile the 
fact that the nonlinear code generates iterates in $[\ell,u]$, while the 
linear solver expects and returns a point in $\real^n$.  The point 
based approximation is formed at $x^k \in [\ell,u]$.  The iterate passed into
the linear solver, $y^k$, is chosen to have the best normal map residual among 
all $y$ such that $\pi_{[\ell,u]}(y) = x^k$.  The point returned from the 
linear solver is termed $x^N$ and is used in a backtracking search instead
of the simple test given as (\ref{descent}).  This search inspects points 
on the following arc parametrized by $t \in (0,1]$:
$$
    \pi_{[\ell,u]}(x^k + t (x^N - x^k)).
$$
Note that, in general, this is not the line segment joining $x^k$ to $x^N$.  
The projection is necessary so that $x^{k+1}$ remains in $[\ell,u]$.
A projected gradient step was only taken if suitable descent did not occur 
for some minimum step length allowed.  Secondly, the gradient of the merit 
function required for (\ref{LSAlg}) was calculated using the formulas detailed
in Propositions \ref{PsiDiff} and \ref{GenJac} and the method described in 
\cite{billups:algorithms}.

Many heuristics have been incorporated into the implementation of PATH 4.x.  
These are documented in the sequel.  Note that previous version of PATH 
\cite{dirkse:robust,dirkse.ferris:path,dirkse.ferris:pathsearch} are
based entirely on the algorithms in \cite{robinson:newtons,ralph:global}
and are globalized by using a pathsearch on the two-norm squared of the 
normal map.  Hence, all of the searches were updating to use the 
Fischer-Burmeister merit function.

\subsection{Crashing Method}
The crashing technique \cite{dirkse.ferris:crash} is used to quickly 
identify an active set from the user-supplied starting point.
At this time, a proximal perturbation scheme 
\cite{billups:algorithms,billups:improving} is used to overcome problems 
with a singular basis matrix.  The proximal perturbation is introduced in
the crash method when the matrix factored is determined to be singular.
The value of the perturbation is based on the current merit function value.

Even if the crash method is turned off, a perturbation is added to the
model if the matrix that would have initially been factored in the crash 
method is numerically singular.  This behavior is extremely useful for 
introducing a perturbation for singular models.

\subsection{Nonmonotone Searches}
The first line of defense against convergence to stationary points is the use 
of a nonmonotone linesearch \cite{grippo.lampariello.ea:nonmonotone,
grippo.lampariello.ea:class, ferris.lucidi:nonmonotone}.  In this case we 
define a reference value, $R^k$ and we use this value in the test for 
sufficient decrease
\[
    \Psi (x^k(t_k)) \leq \real^k - \sigma \nabla \Psi(x^k)^T(x^k - x^k(t_k)).
\]
Depending upon the choice of the reference value, this allows the merit 
function to increase from one iteration to the next.  This strategy can not 
only improve convergence, but can also avoid local minimizers by allowing 
such increases.

We now need to detail the choice of the reference value.  Let
$\{M_1,\ldots,M_m\}$ be a finite set of values initialized to
$\kappa \Psi(x^0)$, where $\kappa$ is used to determine the initial set of 
acceptable merit function values and defaults to 20 in the code.  Note 
$\kappa=1$ indicates that we are not going to allow the merit function to 
increase beyond the initial value.  

Having defined the values of $\{M_1,\ldots,M_m\}$ (where we use $m=10$ by 
default), we can now calculate a reference value.  Assuming that $d^k$ is 
the Newton direction, we define $i_0 = \mbox{argmax }M_i$ and $R^k = M_{i_0}$.
After the nonmonotone linesearch finds $t_k$, we update the memory 
so that $M_{i_0} = \Psi(x^k+t_kd^k)$.  That is, we remove the element 
from the memory having the largest merit function value.

When a gradient step is used, it is beneficial to let $x^k = x^{\mbox{best}}$ 
where $x^{\mbox{best}}$ is the point with the absolute best merit function
value encountered so far.   We then recalculate $d^k = -\nabla \Psi(x^k)$
using the best point and let $R^k = \Psi(x^k)$.  That is, we 
force decrease from the best iterate found whenever a gradient step is 
performed.  After a successful step we set $M_i = \Psi(x^k + t_kd^k)$ for
all $i \in [1, \ldots, m]$.  This prevents future iterates from returning
to the same problem area.

A watchdog strategy \cite{chamberlain.powell.ea:watchdog} is also available 
for use in the code.  The method employed allows steps to be accepted when 
they are ``close'' to the current iterate.  Nonmonotonic decrease is enforced 
every 10 iterations by default.

When the nonmonotone linesearch and projected gradient steps fail, the 
problem is restarted from the initial starting point with a modified set 
of options.  These restarts change the algorithm in the hopes that the 
modified algorithm leads to a solution.  The ordering and nature of the 
restarts were determined by empirical evidence based upon tests performed on 
real-world problems.  

\subsection{Linear Complementarity Subproblems}
All versions of PATH solve a linear complementarity problem each major 
iteration.
Let $M \in \Re^{n\times n}$, $q \in \Re^n$, and $[\ell,u]$ be given.  
$(\bar{z}, \bar{w}, \bar{v})$ solves the linear mixed complementarity problem 
defined by $M$, $q$, and $[\ell,u]$ if and only if it satisfies the following 
constrained system of equations:
\begin{eqnarray}
Mz - w + v + q = 0 \label{eqn:f1}\\
w^T(z - \ell) = 0 \label{eqn:f2} \\
v^T(u - z) = 0 \label{eqn:f3} \\
z \in B, w \in \Re^n_+, v \in \Re^n_+,
\end{eqnarray}
where $x + \infty = \infty$ for all $x \in \Re$ and $0\cdot\infty = 0$
by convention.  A triple, $(\hat{z}, \hat{w}, \hat{v})$, satisfying 
equations \eqnref{f1} - \eqnref{f3} is called a complementary triple.  

The objective of the linear model solver is to construct a path from a given 
complementary triple $(\hat{z}, \hat{w}, \hat{v})$ to a solution 
$(\bar{z}, \bar{w}, \bar{v})$.  The algorithm used to solve the linear problem
is identical to that given in \cite{dirkse:robust}; however, artificial 
variables are incorporated into the model and the residual is scaled by $s$.
The augmented system is then:
\begin{eqnarray}
Mz - w + v + Da + \frac{(1-t)}{s}(sr) + q = 0\\
w^T(z-\ell) = 0 \\
v^T(u-z) = 0 \\
z \in B, w \in \Re^n_+, v \in \Re^n_+, a \equiv 0, t \in [0,1]
\end{eqnarray}
where $r$ is the residual, $t$ is the path parameter, $a$ is a vector
of artificial variables, and $D$ is a permuted diagonal matrix. 

Artificial variables are used to construct an initial invertible basis 
consistent with the given starting point even under
rank deficiency.  The procedure consists of two parts: constructing an
initial guess as to the basis and then recovering from rank deficiency to
obtain an invertible basis.  
The crash technique gives a good approximation to the active
set which is used in the first phase of the algorithm to construct a 
basis by partitioning the variables into three sets:
\begin{enumerate}
\item $W = \{i \in \{1, \ldots, n\} \mid \hat{z}_i = l_i \mbox{ and } \hat{w}_i > 0\}$
\item $V = \{i \in \{1, \ldots, n\} \mid \hat{z}_i = u_i \mbox{ and } \hat{v}_i > 0\}$
\item $Z = \{1, \ldots, n\} \setminus W\cup V$
\end{enumerate}
Since $(\hat{z}, \hat{w}, \hat{v})$ is a complementary triple, 
$Z\cap W\cap V = \emptyset$ and $Z\cup W\cup V = \{1, \ldots, n\}$.
Using the above guess, we can recover an invertible basis consistent with the
starting point by defining $D$ appropriately.  The technique relies upon the 
factorization to tell the linearly dependent rows and columns of the basis 
matrix.  Some of the variables may be nonbasic, but not at their bounds.  
For such variables, the corresponding artificial will be basic.  

We use a modified version of EXPAND \cite{gill.murray.ea:practical*1} 
to perform the ratio test.  Variables are prioritized as follows:
\begin{enumerate}
\item $t$ leaving at its upper bound.
\item Any artificial variable.
\item Any $z$, $w$, or $v$ variable.
\end{enumerate}
If a choice as to the leaving variable can be made while maintaining
numerical stability and sparsity, we choose the variable with the highest 
priority (lowest number above).

When an artificial variable leaves the basis and a $z$-type variable enters, 
we have the choice of either increasing or decreasing that entering
variable because it is nonbasic but not at a bound.  The determination
is made such that $t$ increases and stability is preserved.

If the code is forced to use a ray start at each iteration, then the code 
carries out Lemke's method, which is known 
\cite{cottle.dantzig:complementary} not to cycle.  However, by default, we 
use a regular start to guarantee that the generated path emanates from the 
current iterate.  Under appropriate conditions, this guarantees a decrease 
in the normal map residual.  However, it is then possible for the pivot 
sequence in the linear model to cycle.  To prevent this undesirable
outcome,  we attempt to detect the formation of a cycle with the heuristic that
if a variable enters the basis more that a given number of times, we are
cycling.  The number of times the variable has entered is reset whenever $t$
increases beyond its previous maximum or an artificial variable leaves the
basis.  If cycling is detected, we terminate the linear solver at the largest 
value of $t$ and return this point.

Another heuristic is added when the linear code terminates on a ray.  The
returned point in this case is not the base of the ray.  We move a slight
distance up the ray and return this new point.  The next linear subproblem
should then have a different point-based approximation, which we attempt
to solve.

Since the EXPAND pivot rules are used, some of the variables may be
nonbasic, but slightly infeasible, at the solution.  Whenever the
linear code finishes, the nonbasic variables are put at their bounds 
and the basic variables are recomputed using the current factorization.  
The final step projects the recomputed basic variables onto their bounds.  
This correction procedure helps to find the best possible solution to the 
linear system.

\subsection{Other Features}
Other heuristics are incorporated into the code.  During the first iteration, 
if the linear solver fails to find a Newton point, a Lemke ray start is 
used.  Furthermore, if we repeatedly fail to solve the linear subproblems, an 
advanced ray start will be attempted.  The advanced ray start chooses 
the ``closest'' extreme point of the $[\ell, u]$ polytope and selects a ray 
in the interior of the normal cone at this point.  This technique helps to 
reduce the number of 
pivots required to solve the linear subproblem when compared to a Lemke
ray start.  However, when the basis corresponding to the cell of the
normal manifold selected is not invertible the procedure fails and we resort 
to using a Lemke ray start.  Computational experience has shown this to be an 
effective heuristic and generally results in solving the linear model.  Using 
Lemke ray starts are not the default mode, since typically many more pivots 
are required.

The proximal point perturbation is shrunk each major iteration.  
However, when numerical difficulties are encountered, the value is
increased to a fraction of the current merit function value.

We finally note, that when the merit function fails to show sufficient
decrease over the last 100 iterates, a restart will be performed, as
this tends to indicate that we are close to a stationary point.

\section{Computational Results}\label{sec:path-results}
In the following tables, we report the number of successes and failures 
(the latter in parenthesis) of the new code, PATH 4.x, compared to
the PATH 2.9 code \cite{dirkse:robust} from all starting points in the 
GAMSLIB, MCPLIB \cite{dirkse.ferris:mcplib}, and NETLIB \cite{gay:electronic}
collections of test problems.  The NETLIB problems are originally linear 
programs. We formed the optimality conditions for these linear programs 
and solved the resulting complementarity problems for the NETLIB tests.

\begin{table}
\centering
\caption{Comparison of PATH 4.x and PATH 2.9 on GAMSLIB}\label{tab:path-gamslib}
\begin{tabular}{l|c|cc}
Problem & & Without Restarts & With Restarts \\
        & PATH 2.9 & PATH 4.x & PATH 4.x \\
\hline
cafemge & 2 (0) &	1 (1) &	2 (0) \\
cammcp & 1 (0) &	1 (0) &	1 (0) \\
cammge & 1 (0) &	1 (0) &	1 (0) \\
cirimge & 6 (0) &	6 (0) &	6 (0) \\
co2mge & 6 (0) &	6 (0) &	6 (0) \\
dmcmge & 2 (0) &	2 (0) &	2 (0) \\
ers82mcp & 1 (0) &	1 (0) &	1 (0) \\
etamge & 1 (0) &	1 (0) &	1 (0) \\
finmge & 4 (0) &	4 (0) &	4 (0) \\
gemmcp & 4 (0) &	4 (0) &	4 (0) \\
gemmge & 4 (0) &	4 (0) &	4 (0) \\
hansmcp & 1 (0) &	1 (0) &	1 (0) \\
hansmge & 1 (0) &	1 (0) &	1 (0) \\
harkmcp & 4 (0) &	4 (0) &	4 (0) \\
harmge & 6 (0) &	6 (0) &	6 (0) \\
kehomge & 3 (0) &	3 (0) &	3 (0) \\
kormcp & 1 (0) &	1 (0) &	1 (0) \\
mr5mcp & 1 (0) &	1 (0) &	1 (0) \\
nsmge & 1 (0) &	1 (0) &	1 (0) \\
oligomcp & 1 (0) &	1 (0) &	1 (0) \\
sammge & 16 (0) &	16 (0) &	16 (0) \\
scarfmcp & 1 (0) &	1 (0) &	1 (0) \\
scarfmge & 4 (0) &	4 (0) &	4 (0) \\
shovmge & 4 (0) &	4 (0) &	4 (0) \\
threemge & 9 (0) &	9 (0) &	9 (0) \\
transmcp & 3 (0) &	3 (0) &	3 (0) \\
two3mcp & 2 (0) &	2 (0) &	2 (0) \\
unstmge & 1 (0) &	1 (0) &	1 (0) \\
vonthmcp & 1 (0) &	0 (1) &	1 (0) \\
vonthmge & 1 (0) &	1 (0) &	1 (0) \\
wallmcp & 1 (0) &	1 (0) &	1 (0) \\
\hline
Total & 94 (0) &	92 (2) &	94 (0)
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Comparison of PATH 4.x and PATH 2.9 on MCPLIB}\label{tab:path-mcplib}
\begin{tabular}{l|c|cc}
Problem & & Without Restarts & With Restarts \\
        & PATH 2.9 & PATH 4.x & PATH 4.x \\
\hline
asean9a & 1 (0) &	1 (0) &	1 (0) \\
badfree & 0 (1) &	1 (0) &	1 (0) \\
bai\_haung & 1 (0) &	1 (0) &	1 (0) \\
bert\_oc & 4 (0) &	4 (0) &	4 (0) \\
bertsekas & 6 (0) &	6 (0) &	6 (0) \\
billups & 0 (3) &	0 (3) &	0 (3) \\
bishop & 1 (0) &	1 (0) &	1 (0) \\
box & 352 (9) &	361 (0) &	361 (0) \\
bratu & 1 (0) &	1 (0) &	1 (0) \\
cammcf & 0 (1) &	0 (1) &	0 (1) \\
carbon & 6 (1) &	6 (1) &	6 (1) \\
cgereg & 19 (0) &	19 (0) &	19 (0) \\
choi & 1 (0) &	1 (0) &	1 (0) \\
colvdual & 4 (0) &	2 (2) &	4 (0) \\
colvnlp & 6 (0) &	6 (0) &	6 (0) \\
congest & 1 (0) &	1 (0) &	1 (0) \\
cycle & 1 (0) &	1 (0) &	1 (0) \\
degen & 1 (0) &	1 (0) &	1 (0) \\
denmark & 34 (4) &	38 (0) &	38 (0) \\
dirkse1 & 0 (1) &	1 (1) &	1 (1) \\
duopoly & 0 (1) &	0 (1) &	0 (1) \\
eckstein & 0 (1) &	1 (0) &	1 (0) \\
ehl\_k40 & 2 (1) &	2 (1) &	3 (0) \\
ehl\_k60 & 3 (0) &	2 (1) &	3 (0) \\
ehl\_k80 & 3 (0) &	3 (0) &	3 (0) \\
ehl\_kost & 3 (0) &	3 (0) &	3 (0) \\
electric & 1 (0) &	0 (1) &	1 (0) \\
endog & 0 (3) &	23 (0) &	23 (0) \\
eppa & 8 (0) &	8 (0) &	8 (0) \\
eta2100 & 1 (0) &	1 (0) &	1 (0) \\
exemptions & 181 (5) &	186 (0) &	186 (0) \\
explcp & 1 (0) &	1 (0) &	1 (0) \\
exros & 5 (0) &	2 (3) &	5 (0)
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Comparison of PATH 4.x and PATH 2.9 on MCPLIB (cont.)}\label{tab:path-mcplib2}
\begin{tabular}{l|c|cc}
Problem & & Without Restarts & With Restarts \\
        & PATH 2.9 & PATH 4.x & PATH 4.x \\
\hline
finance1 & 21 (0) &	21 (0) &	21 (0) \\
finance2 & 10 (0) &	10 (0) &	10 (0) \\
finance3 & 20 (0) &	20 (0) &	20 (0) \\
fixedpt & 0 (1) &	0 (1) &	1 (0) \\
forcebsm & 1 (0) &	1 (0) &	1 (0) \\
forcedsa & 1 (0) &	1 (0) &	1 (0) \\
freebert & 7 (0) &	7 (0) &	7 (0) \\
fried7 & 5 (0) &	5 (0) &	5 (0) \\
fried8 & 4 (1) &	2 (3) &	3 (2) \\
gafni & 3 (0) &	3 (0) &	3 (0) \\
games & 22 (3) &	23 (2) &	25 (0) \\
gei & 2 (0) &	2 (0) &	2 (0) \\
golanmcp & 1 (0) &	1 (0) &	1 (0) \\
hanskoop & 10 (0) &	10 (0) &	10 (0) \\
hansmcf & 1 (0) &	1 (0) &	1 (0) \\
hanson & 1 (0) &	1 (0) &	1 (0) \\
ho & 5 (0) &	5 (0) &	5 (0) \\
hydroc06 & 1 (0) &	1 (0) &	1 (0) \\
hydroc20 & 1 (0) &	1 (0) &	1 (0) \\
jel & 2 (0) &	2 (0) &	2 (0) \\
jmu & 1 (0) &	0 (1) &	1 (0) \\
josephy & 8 (0) &	8 (0) &	8 (0) \\
keyzer & 4 (2) &	5 (1) &	6 (0) \\
kojshin & 8 (0) &	8 (0) &	8 (0) \\
kyh-scale & 1 (1) &	0 (2) &	1 (1) \\
kyh & 0 (2) &	2 (0) &	2 (0) \\
leyffer & 1 (0) &	1 (0) &	1 (0) \\
lincont & 1 (0) &	1 (0) &	1 (0) \\
lstest & 0 (1) &	0 (1) &	0 (1) \\
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Comparison of PATH 4.x and PATH 2.9 on MCPLIB (cont.)}\label{tab:path-mcplib3}
\begin{tabular}{l|c|cc}
Problem & & Without Restarts & With Restarts \\
        & PATH 2.9 & PATH 4.x & PATH 4.x \\
\hline
markusen & 18 (0) &	18 (0) &	18 (0) \\
mathinum & 6 (0) &	6 (0) &	6 (0) \\
mathisum & 6 (0) &	6 (0) &	6 (0) \\
methan08 & 1 (0) &	1 (0) &	1 (0) \\
mr5mcf & 1 (0) &	1 (0) &	1 (0) \\
mrt & 1 (0) &	1 (0) &	1 (0) \\
multi-v1 & 1 (0) &	1 (0) &	1 (0) \\
multi-v2 & 1 (0) &	1 (0) &	1 (0) \\
multi-v3 & 1 (0) &	1 (0) &	1 (0) \\
munson3 & 1 (0) &	0 (1) &	1 (0) \\
munson4 & 3 (0) &	3 (0) &	3 (0) \\
nash & 4 (0) &	4 (0) &	4 (0) \\
ne-hard & 0 (1) &	1 (0) &	1 (0) \\
nepal1 & 2 (1) &	0 (3) &	2 (1) \\
nepal2 & 0 (1) &	1 (0) &	1 (0) \\
obstacle & 8 (0) &	8 (0) &	8 (0) \\
olg & 1 (0) &	1 (0) &	1 (0) \\
opening & 0 (2) &	2 (0) &	2 (0) \\
opt\_cont & 1 (0) &	1 (0) &	1 (0) \\
opt\_cont127 & 1 (0) &	1 (0) &	1 (0) \\
opt\_cont255 & 1 (0) &	1 (0) &	1 (0) \\
opt\_cont31 & 1 (0) &	1 (0) &	1 (0) \\
opt\_cont511 & 1 (0) &	1 (0) &	1 (0) \\
pgvon105 & 3 (3) &	5 (1) &	5 (1) \\
pgvon106 & 2 (4) &	5 (1) &	5 (1) \\
pies & 1 (0) &	1 (0) &	1 (0) \\
pizer & 3 (0) &	3 (0) &	3 (0) \\
powell & 5 (1) &	6 (0) &	6 (0) \\
powell\_mcp & 6 (0) &	6 (0) &	6 (0) \\
qp & 1 (0) &	1 (0) &	1 (0) \\
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Comparison of PATH 4.x and PATH 2.9 on MCPLIB (cont.)}\label{tab:path-mcplib4}
\begin{tabular}{l|c|cc}
Problem & & Without Restarts & With Restarts \\
        & PATH 2.9 & PATH 4.x & PATH 4.x \\
\hline
ralph & 6 (1) &	7 (0) &	7 (0) \\
renger & 0 (1) &	0 (1) &	0 (1) \\
romer & 0 (1) &	1 (1) &	2 (0) \\
romer2 & 1 (6) &	0 (7) &	1 (6) \\
romer3 & 0 (1) &	0 (1) &	0 (1) \\
runge & 7 (0) &	7 (0) &	7 (0) \\
scarfanum & 4 (0) &	4 (0) &	4 (0) \\
scarfasum & 4 (0) &	4 (0) &	4 (0) \\
scarfbnum & 2 (0) &	2 (0) &	2 (0) \\
scarfbsum & 1 (1) &	2 (0) &	2 (0) \\
shubik & 43 (4) &	43 (4) &	47 (0) \\
simple-ex & 0 (1) &	0 (1) &	0 (1) \\
simple-red & 1 (0) &	1 (0) &	1 (0) \\
spill1 & 0 (0) &	1 (0) &	1 (0) \\
spill2 & 1 (0) &	0 (1) &	1 (0) \\
spill3 & 0 (1) &	0 (1) &	1 (0) \\
spill4 & 0 (1) &	1 (0) &	1 (0) \\
spill5 & 0 (1) &	1 (0) &	1 (0) \\
sppe & 3 (0) &	3 (0) &	3 (0) \\
tinloi & 58 (6) &	56 (8) &	64 (0) \\
tinloilp & 4 (0) &	4 (0) &	4 (0) \\
tinsmall & 63 (1) &	63 (1) &	64 (0) \\
titan & 1 (1) &	1 (1) &	1 (1) \\
tobin & 4 (0) &	4 (0) &	4 (0) \\
tqbilat & 1 (0) &	1 (0) &	1 (0) \\
trade12 & 2 (0) &	2 (0) &	2 (0) \\
trafelas & 2 (0) &	2 (0) &	2 (0) \\
trig & 3 (0) &	3 (0) &	3 (0) \\
uruguay & 4 (0) &	4 (0) &	4 (0) \\
venables & 1 (1) &	2 (0) &	2 (0) \\
vonthmcf & 1 (0) &	1 (0) &	1 (0) \\
water & 2 (1) &	2 (1) &	3 (0) \\
xu1 & 7 (0) &	7 (0) &	7 (0) \\
xu2 & 7 (0) &	7 (0) &	7 (0) \\
xu3 & 7 (0) &	7 (0) &	7 (0) \\
xu4 & 7 (0) &	7 (0) &	7 (0) \\
xu5 & 7 (0) &	7 (0) &	7 (0) \\
\hline
Total & 1124 (85) &	1171 (60) &	1207 (24)
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Comparison of PATH 4.x and PATH 2.9 on NETLIB}\label{tab:path-netlib}
\begin{tabular}{l|c|cc}
Problem & & Without Restarts & With Restarts \\
        & PATH 2.9 & PATH 4.x & PATH 4.x \\
\hline
adlittle & 1 (0) &	1 (0) &	1 (0) \\
afiro & 1 (0) &	1 (0) &	1 (0) \\
agg & 1 (0) &	1 (0) &	1 (0) \\
agg2 & 1 (0) &	1 (0) &	1 (0) \\
agg3 & 1 (0) &	1 (0) &	1 (0) \\
bandm & 1 (0) &	1 (0) &	1 (0) \\
beaconfd & 1 (0) &	1 (0) &	1 (0) \\
blend & 1 (0) &	1 (0) &	1 (0) \\
bnl1 & 1 (0) &	1 (0) &	1 (0) \\
bnl2 & 0 (1) &	1 (0) &	1 (0) \\
boeing1 & 1 (0) &	1 (0) &	1 (0) \\
boeing2 & 1 (0) &	1 (0) &	1 (0) \\
bore3d & 1 (0) &	1 (0) &	1 (0) \\
brandy & 1 (0) &	1 (0) &	1 (0) \\
capri & 1 (0) &	1 (0) &	1 (0) \\
cycle & 0 (1) &	0 (1) &	0 (1) \\
czprob & 1 (0) &	1 (0) &	1 (0) \\
d2q06c & 0 (1) &	0 (1) &	0 (1) \\
d6cube & 1 (0) &	1 (0) &	1 (0) \\
degen2 & 1 (0) &	1 (0) &	1 (0) \\
degen3 & 1 (0) &	1 (0) &	1 (0) \\
dfl001 & 0 (1) &	0 (1) &	0 (1) \\
e226 & 1 (0) &	1 (0) &	1 (0) \\
etamacro & 1 (0) &	1 (0) &	1 (0) \\
fffff800 & 0 (1) &	1 (0) &	1 (0) \\
finnis & 1 (0) &	1 (0) &	1 (0) \\
fit1d & 1 (0) &	1 (0) &	1 (0) \\
fit1p & 1 (0) &	1 (0) &	1 (0) \\
fit2d & 0 (1) &	0 (1) &	0 (1) \\
fit2p & 0 (1) &	0 (1) &	0 (1) \\
forplan & 1 (0) &	1 (0) &	1 (0) \\
ganges & 1 (0) &	1 (0) &	1 (0) \\
gfrd-pnc & 0 (1) &	1 (0) &	1 (0) \\
greenbea & 0 (1) &	0 (1) &	0 (1) \\
greenbeb & 0 (1) &	0 (1) &	0 (1) \\
grow15 & 1 (0) &	1 (0) &	1 (0) \\
grow22 & 1 (0) &	1 (0) &	1 (0) \\
grow7 & 1 (0) &	1 (0) &	1 (0) \\
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Comparison of PATH 4.x and PATH 2.9 on NETLIB (cont.)}\label{tab:path-netlib2}
\begin{tabular}{l|c|cc}
Problem & & Without Restarts & With Restarts \\
        & PATH 2.9 & PATH 4.x & PATH 4.x \\
\hline
israel & 1 (0) &	1 (0) &	1 (0) \\
kb2 & 1 (0) &	1 (0) &	1 (0) \\
lotfi & 1 (0) &	1 (0) &	1 (0) \\
nesm & 0 (1) &	0 (1) &	0 (1) \\
perold & 0 (1) &	1 (0) &	1 (0) \\
pilot4 & 0 (1) &	0 (1) &	0 (1) \\
sc105 & 1 (0) &	1 (0) &	1 (0) \\
sc205 & 1 (0) &	1 (0) &	1 (0) \\
sc50a & 1 (0) &	1 (0) &	1 (0) \\
sc50b & 1 (0) &	1 (0) &	1 (0) \\
scagr25 & 1 (0) &	1 (0) &	1 (0) \\
scagr7 & 1 (0) &	1 (0) &	1 (0) \\
scfxm1 & 1 (0) &	1 (0) &	1 (0) \\
scfxm2 & 1 (0) &	1 (0) &	1 (0) \\
scfxm3 & 1 (0) &	1 (0) &	1 (0) \\
scorpion & 1 (0) &	1 (0) &	1 (0) \\
scrs8 & 1 (0) &	1 (0) &	1 (0) \\
scsd1 & 1 (0) &	1 (0) &	1 (0) \\
scsd6 & 1 (0) &	1 (0) &	1 (0) \\
scsd8 & 1 (0) &	1 (0) &	1 (0) \\
sctap1 & 1 (0) &	1 (0) &	1 (0) \\
sctap2 & 1 (0) &	1 (0) &	1 (0) \\
sctap3 & 1 (0) &	1 (0) &	1 (0) \\
share1b & 1 (0) &	1 (0) &	1 (0) \\
share2b & 1 (0) &	1 (0) &	1 (0) \\
shell & 1 (0) &	1 (0) &	1 (0) \\
ship04l & 1 (0) &	1 (0) &	1 (0) \\
ship04s & 1 (0) &	1 (0) &	1 (0) \\
ship08l & 1 (0) &	1 (0) &	1 (0) \\
ship08s & 1 (0) &	1 (0) &	1 (0) \\
ship12l & 0 (1) &	1 (0) &	1 (0) \\
ship12s & 1 (0) &	1 (0) &	1 (0) \\
stair & 1 (0) &	1 (0) &	1 (0) \\
stocfor1 & 1 (0) &	1 (0) &	1 (0) \\
stocfor2 & 1 (0) &	1 (0) &	1 (0) \\
tuff & 0 (1) &	1 (0) &	1 (0) \\
wood1p & 0 (1) &	1 (0) &	1 (0) \\
\hline
Total & 59 (16) &	66 (9) &	66 (9) \\
\end{tabular}
\end{table}

In the reported results, PATH 4.x code is able to carry out at most 5 
projected gradient steps when the direction provided by the linear 
subproblem is poor.  Since the merit function in PATH 2.9 is nonsmooth, 
there is no possibility of carrying out a similar scheme in this code.  
The projected gradient steps enable progress to be made in PATH 4.x, which 
cannot occur in PATH 2.9.  We note that a total of 361 projected gradient 
steps were taken during the course of 1400 runs, indicating the preference to 
take such steps only as a last resort (i.e. only after the nonmonotone search 
and the watchdog strategy fail).  Several heuristic procedures were described 
and tested in \cite{dirkse.ferris:crash}; all of these appear not to be 
beneficial to PATH 4.x. 

PATH 2.9 did not encounter all of the subproblems for \texttt{endog}, 
\texttt{romer}, and \texttt{spill1} because of failures.  Thus, there
is a discrepancy between the totals reported in the MCPLIB test set.
However, we can see from the results that PATH 4.x is considerably more
reliable than PATH 2.9 even when restarts are not allowed.  Note that there 
is  no way to turn off the ad-hoc heuristics contained in PATH 2.9, so the
appropriate comparison is with the version of PATH 4.x allowing restarts.

To conclude this section, \tabref{path-complete} and 
\tabref{path-complete2} provide more complete information for PATH 4.x on a 
subset of the test problems.  These problems were selected by taking every 
10th run in an alphabetical ordering of all the test problems considered.  
However, we only report the first three instances of any given problem.  We 
believe this is an unbiased sample of the test results.
 
\begin{table}
\centering\small
\caption{Selected full results for PATH 4.x}\label{tab:path-complete}
\begin{tabular}{llc|ccccc}
Library & Problem & SP & MI & CI & R & PG & T \\
\hline
gamslib & cirimge & 6 & 4 & 0 & 0 & 0 & 0.0 (0.0) \\
gamslib & etamge & 1 & 7 & 0 & 0 & 0 & 0.1 (0.1) \\
gamslib & gemmge & 3 & 5 & 0 & 0 & 0 & 0.3 (0.3) \\
gamslib & harmge & 2 & 4 & 0 & 0 & 0 & 0.0 (0.0) \\
gamslib & nsmge & 1 & 9 & 0 & 0 & 0 & 0.2 (0.1) \\
gamslib & sammge & 11 & 3 & 0 & 0 & 0 & 0.0 (0.0) \\
gamslib & scarfmge & 2 & 6 & 0 & 0 & 0 & 0.0 (0.0) \\
gamslib & threemge & 6 & 4 & 0 & 0 & 0 & 0.0 (0.0) \\
gamslib & two3mcp & 2 & 3 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & bert\_oc & 3 & 0 & 3 & 0 & 0 & 0.6 (0.3) \\
mcplib & billups & 3 & 16 & 2 & 3 & 0 & * (*) \\
mcplib & box & 9 & 3 & 0 & 0 & 0 & 0.0 (0.0) \\
mcplib & box & 19 & 4 & 0 & 0 & 0 & 0.0 (0.0) \\
mcplib & box & 29 & 3 & 0 & 0 & 0 & 0.0 (0.0) \\
mcplib & carbon & 7 & 3 & 0 & 0 & 0 & 0.4 (0.4) \\
mcplib & cgereg & 12 & 1 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & cgereg & 22 & 12 & 1 & 0 & 0 & 0.2 (0.2) \\
mcplib & colvnlp & 5 & 6 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & denmark & 9 & 10 & 0 & 0 & 0 & 9.0 (8.5) \\
mcplib & denmark & 26 & 12 & 0 & 0 & 0 & 4.2 (2.9) \\
mcplib & denmark & 42 & 8 & 0 & 0 & 0 & 7.8 (10.6) \\
mcplib & ehl\_k60 & 1 & 4 & 1 & 0 & 0 & 0.3 (0.2) \\
mcplib & endog & 2 & 14 & 0 & 0 & 0 & 26.80 (*) \\
mcplib & endog & 12 & 2 & 0 & 0 & 0 & 7.18 (*) \\
mcplib & endog & 22 & 1 & 0 & 0 & 0 & 5.37 (*) \\
mcplib & eppa & 8 & 4 & 0 & 0 & 0 & 1.7 (1.6) \\
mcplib & exemptions & 10 & 7 & 0 & 0 & 0 & 3.6 (1.3) \\
mcplib & exemptions & 20 & 4 & 0 & 0 & 0 & 2.3 (1.6) \\
mcplib & exemptions & 30 & 8 & 0 & 0 & 0 & 3.8 (2.1) \\
mcplib & exros & 2 & 115 & 1 & 1 & 0 & 14.2 (0.2) \\
mcplib & finance1 & 7 & 0 & 2 & 0 & 0 & 0.0 (0.0) \\
mcplib & finance1 & 17 & 0 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & finance2 & 6 & 1 & 2 & 0 & 0 & 0.1 (0.0) \\
mcplib & finance3 & 6 & 0 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & finance3 & 16 & 1 & 1 & 0 & 0 & 0.1 (0.0) \\
mcplib & freebert & 3 & 3 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & fried8 & 1 & 191 & 1 & 2 & 7 & 3.86 (*) \\
mcplib & games & 3 & 42 & 1 & 0 & 0 & 0.2 (0.1) \\
mcplib & games & 13 & 7 & 0 & 0 & 0 & 0.1 (0.0) \\
mcplib & games & 23 & 5 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & hanskoop & 5 & 5 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & ho & 3 & 2 & 0 & 0 & 0 & 0.0 (0.0) \\
\end{tabular}
\end{table}

\begin{table}
\centering\small
\caption{Selected full results for PATH 4.x (cont.)}\label{tab:path-complete2}
\begin{tabular}{llc|ccccc}
Library & Problem & SP & MI & CI & R & PG & T \\
\hline
mcplib & josephy & 3 & 14 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & keyzer & 5 & 6 & 0 & 0 & 0 & 0.0 (0.0) \\
mcplib & kyh-scale & 1 & 202 & 7 & 3 & 13 & * (0.2) \\
mcplib & markusen & 6 & 3 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & markusen & 24 & 5 & 0 & 0 & 0 & 0.0 (0.0) \\
mcplib & mathinum & 6 & 6 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & multi-v1 & 1 & 15 & 0 & 0 & 0 & 0.1 (0.0) \\
mcplib & nash & 4 & 2 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & obstacle & 5 & 1 & 1 & 0 & 0 & 9.4 (1.0) \\
mcplib & opt\_cont31 & 1 & 1 & 4 & 0 & 0 & 0.4 (0.2) \\
mcplib & pgvon106 & 3 & 18 & 1 & 0 & 0 & 0.3 (0.2) \\
mcplib & powell & 3 & 6 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & qp & 1 & 1 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & romer & 2 & 25 & 0 & 0 & 1 & 1.33 (*) \\
mcplib & runge & 2 & 23 & 0 & 0 & 0 & 0.0 (0.0) \\
mcplib & scarfasum & 1 & 4 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & shubik & 3 & 9 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & shubik & 14 & 13 & 2 & 0 & 0 & 0.0 (*) \\
mcplib & shubik & 24 & 2 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & spill4 & 1 & 29 & 0 & 0 & 2 & 25.96 (*) \\
mcplib & tinloi & 6 & 1 & 2 & 0 & 0 & 0.1 (0.0) \\
mcplib & tinloi & 16 & 137 & 1 & 1 & 0 & 2.6 (*) \\
mcplib & tinloi & 26 & 1 & 3 & 0 & 0 & 0.1 (0.0) \\
mcplib & tinloilp & 2 & 3 & 2 & 0 & 0 & 114.1 (141.2) \\
mcplib & tinsmall & 8 & 1 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & tinsmall & 18 & 2 & 1 & 0 & 0 & 0.0 (0.1) \\
mcplib & tinsmall & 28 & 2 & 1 & 0 & 0 & 0.0 (0.1) \\
mcplib & tobin & 2 & 7 & 2 & 0 & 0 & 0.0 (0.0) \\
mcplib & trig & 3 & 68 & 1 & 0 & 1 & 0.1 (0.2) \\
mcplib & water & 5 & 16 & 0 & 0 & 0 & 3.3 (30.0) \\
mcplib & xu2 & 3 & 4 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & xu3 & 6 & 4 & 1 & 0 & 0 & 0.0 (0.0) \\
mcplib & xu5 & 2 & 5 & 1 & 0 & 0 & 0.0 (0.0) \\
netlib & agg3 & 1 & 22 & 2 & 0 & 0 & 1.1 (1.1) \\
netlib & capri & 1 & 16 & 17 & 0 & 0 & 1.6 (2.6) \\
netlib & fffff800 & 1 & 27 & 0 & 0 & 0 & 15.8 (*) \\
netlib & greenbeb & 1 & 9 & 29 & 0 & 0 & * (*) \\
netlib & sc105 & 1 & 9 & 1 & 0 & 0 & 0.1 (0.1) \\
netlib & scrs8 & 1 & 15 & 1 & 0 & 0 & 3.8 (4.0) \\
netlib & ship04l & 1 & 13 & 0 & 0 & 0 & 2.9 (2.6) \\
netlib & wood1p & 1 & 4 & 0 & 0 & 0 & 465.18 (*) \\
\end{tabular}
\end{table}

The columns of these tables indicate the starting point number 
(SP), number of major iterations (MI), crash iterations (CI), restarts (R), 
and projected gradient steps (PG) taken.  The final column of this table 
reports the time for PATH 4.x in seconds, with the time for PATH 2.9 added in 
parentheses.  All runs were carried out on a Sun UltraSparc 330 MHz processor 
with 768 MB RAM.  A ``*'' indicates failure of the method.  Both codes were 
compiled with the same compiler and options.

It is hard to draw firm conclusions from \tabref{path-complete} and
\tabref{path-complete2}.  The solution times of both algorithms are 
comparable, with some smaller times for PATH 4.x and some for PATH 2.9.  
There are only 6 problems in this subset which use projected gradient steps 
and restarts, the vast majority of them being solved without invoking these 
strategies.  Overall, the theoretical extensions outlined and implemented 
result in the improved robustness of PATH 4.x over PATH 2.9 without any 
substantial changes in accuracy or speed.

\section{User Documentation}\label{sec:user}
The previous sections discussed the theoretical foundations of PATH 4.x,
presented the implementation features, and demonstrated the reliability
of the code.  We now detail the PATH 4.x code from the perspective of a user.  
The log file informs a user about any progress made by the code.  If the 
algorithm fails to solve a particular model, the log file can be used to 
track the reason for the failure.  We document the output typically produced 
by PATH 4.x in \secref{anatomy}, which assumes that the modeler is attempting
to solve the \texttt{transport} model developed in \secref{applications}.
\secref{path-options} discusses the available options used to modify the
behavior of the algorithm.  These options can be set by a modeler to 
improve the performance of the code for their particular complementarity
problem.

\subsection{Anatomy of the Log File}\label{sec:anatomy}
We now describe the behavior of the PATH 4.x algorithm in terms of
the output produced for the \texttt{transport} model.  
An example of the log for a particular run is given in \figref{logfile} 
and \figref{logfile2}.

\begin{figure}
\begin{smallverb}
Path 4.4a (Sat Feb 26 09:38:08 2000)

INITIAL POINT STATISTICS
Maximum of X. . . . . . . . . . -0.0000e+00 var: (x.seattle.new-york)
Maximum of F. . . . . . . . . .  6.0000e+02 eqn: (supply.san-diego)
Maximum of Grad F . . . . . . .  1.0000e+00 eqn: (demand.new-york)
                                            var: (x.seattle.new-york)

INITIAL JACOBIAN NORM STATISTICS
Maximum Row Norm. . . . . . . .  3.0000e+00 eqn: (supply.seattle)
Minimum Row Norm. . . . . . . .  2.0000e+00 eqn: (rational.seattle.new-york)
Maximum Column Norm . . . . . .  3.0000e+00 var: (p_supply.seattle)
Minimum Column Norm . . . . . .  2.0000e+00 var: (x.seattle.new-york)

Crash Log
major  func  diff  size  residual    step       prox   (label)
    0     0             1.0416e+03             0.0e+00 (demand.new-york)
    1     1     3     3 1.0029e+03  1.0e+00    1.0e+01 (demand.new-york)
pn_search terminated: no basis change.
\end{smallverb}
\caption{Log File from PATH for solving \texttt{transmcp}}
\label{fig:logfile}
\end{figure}

\begin{figure}
\begin{smallverb}
Major Iteration Log
major minor  func  grad  residual    step  type prox    inorm  (label)
    0     0     2     2 1.0029e+03           I 9.0e+00 6.2e+02 (demand.new-york)
    1     1     3     3 8.3096e+02  1.0e+00 SO 3.6e+00 4.5e+02 (demand.new-york)

...

   15     2    17    17 1.3972e-09  1.0e+00 SO 4.8e-06 1.3e-09 (demand.chicago)

FINAL STATISTICS
Inf-Norm of Complementarity . .  1.4607e-08 eqn: (rational.seattle.chicago)
Inf-Norm of Normal Map. . . . .  1.3247e-09 eqn: (demand.chicago)
Inf-Norm of Minimum Map . . . .  1.3247e-09 eqn: (demand.chicago)
Inf-Norm of Fischer Function. .  1.3247e-09 eqn: (demand.chicago)
Inf-Norm of Grad Fischer Fcn. .  1.3247e-09 eqn: (rational.seattle.chicago)

FINAL POINT STATISTICS
Maximum of X. . . . . . . . . .  3.0000e+02 var: (x.seattle.chicago)
Maximum of F. . . . . . . . . .  5.0000e+01 eqn: (supply.san-diego)
Maximum of Grad F . . . . . . .  1.0000e+00 eqn: (demand.new-york)
                                            var: (x.seattle.new-york)

 ** EXIT - solution found.

Major Iterations. . . . 15
Minor Iterations. . . . 31
Restarts. . . . . . . . 0
Crash Iterations. . . . 1
Gradient Steps. . . . . 0
Function Evaluations. . 17
Gradient Evaluations. . 17
Total Time. . . . . . . 0.020000
Residual. . . . . . . . 1.397183e-09
\end{smallverb}
\caption{Log File from PATH for solving \texttt{transmcp} (cont.)}
\label{fig:logfile2}
\end{figure}

After some basic memory allocation and problem checking, PATH 4.x checks if the 
modeler required an option file to be read.  In this particular example, we 
did not request that an option file be used.  However, if PATH 4.x is directed 
to read an option file, then the following output is generated after the 
PATH 4.x banner.
\single\begin{verbatim}
     Reading options file PATH.OPT
      > output_linear_model yes;
     Options: Read: Line 2 invalid: hi_there;
     Read of options file complete.
\end{verbatim}
\double
If the option reader encounters an invalid option (as above), it
reports this but carries on executing the algorithm.  
The user specifies option name, value pairs in the file.  Note that
only the first three letters of each component of the option name is
relevant.  Therefore, the above option could be set with the
string \texttt{out\_lin\_mod yes}.  A complete listing of the
options available are given in \secref{path-options}.  Following the 
setting of the options, the algorithm starts working on the model.

\subsubsection{Preprocessing}
The first phase of PATH 4.x is to preprocess the model to
reduce the size and complexity.  However, the transportation model given 
to PATH 4.x is not reducible.  For models that are reducible, the preprocessor 
reports its finding with output to the log file.  An example of the 
preprocessing for the \texttt{forcebsm} model is presented below.
\single\begin{verbatim}
     Zero:     0 Single:   112 Double:     0 Forced:     0
     Preprocessed size: 72
\end{verbatim}
\double
The keywords indicate the type of elimination performed, while the
corresponding integer tells the number of constraints where the rule
was used.  An indication of the size of the preprocessed model is
also reported.  See \chapref{preprocessing} for a complete description
of the preprocessor.

On exit from the algorithm, we generate a solution for the original problem
during a postsolve step.  Following the postsolve, the residual of the 
solution generated for the original model is reported. 
\single\begin{verbatim}
     Postsolved residual: 1.0518e-10
\end{verbatim}
\double
This number should be approximately the same as the final residual reported 
for the presolved model.

\subsubsection{Diagnostic Information}
Following the preprocessing step, diagnostic information is generated 
at the starting point.  This information is contained under the
heading ``INITIAL POINT STATISTICS''.  Included is information about the 
initial point and function evaluation.  The log file tells the value of the 
largest element of the starting point and the variable where it occurs.  
Similarly, the maximum function value is displayed along with the name
of the equation producing it.  The maximum element in the gradient is also 
presented with the equation and variable where it is located.

The second block, ``INITIAL JACOBIAN NORM STATISTICS'', provides more 
information about the Jacobian at the starting point.  This information can 
be used to help scale the model as detailed in \chapref{diagnostic}.

\subsubsection{Crash Log}
The main computation in PATH 4.x begins with the crash procedure which 
attempts to quickly determine which of the inequalities should be active.
This procedure is documented fully in \cite{dirkse.ferris:crash}, and
an example of the Crash Log can be seen in \figref{logfile}.
The first column of the crash log is just a label indicating the current
iteration number, the second gives an indication of how many function
evaluations have been performed so far.  Note that precisely one Jacobian 
(gradient) evaluation is performed per crash iteration. The number of changes 
to the active set between iterations of the crash procedure is shown under
the ``diff'' column.  The crash procedure terminates if this becomes
small.  Each iteration of this procedure involves a factorization of a matrix
whose size is shown in the next column.  The residual is a measure of
how far the current iterate is from satisfying the complementarity
conditions; it is zero at a solution.  The column ``step'' corresponds 
to the step length taken in this iteration - ideally this should be 1.
If the factorization fails, then the matrix is perturbed by an identity
matrix scaled by the value indicated in the ``prox'' column.  
The ``label'' column indicates which row in the model is furthest away
from satisfying the complementarity conditions.  Typically, relatively few 
crash 
iterations are performed.  \secref{path-options} gives the options used
to change the behavior of these steps.

\subsubsection{Major Iteration Log}
After the crash is completed, the main algorithm starts as indicated by the
``Major Iteration Log'' flag (see \figref{logfile2}).  The columns that have 
the same labels as in the crash log have precisely the same meaning described 
above.  However, there are some new columns that we now explain.  Each major 
iteration attempts to solve a linear mixed complementarity problem using a 
pivotal method that is a generalization of Lemke's method 
\cite{lemke.howson:equilibrium}.  The number of pivots performed \texttt{per} 
major iteration is given in the ``minor'' column.  

The ``grad'' column gives the cumulative number of Jacobian
evaluations used; typically one evaluation is performed per iteration.
The ``inorm'' column gives the value of the error in satisfying the equation 
indicated in the ``label'' column.

At each iteration of the algorithm, several different step types can be taken,
due to the use of nonmonotone searches 
\cite{dirkse.ferris:pathsearch,ferris.lucidi:nonmonotone}, which are used to
improve robustness.  In order to help the PATH 4.x user, we have added two code 
letters indicating the return code from the linear solver and the step type 
to the log file.  \tabref{7.5} explains the return codes for the linear solver
and \tabref{8} explains the meaning of each step type.  The ideal output in 
this column is ``SO'', with ``SD'' and ``SB'' also being reasonable.  Codes 
different from these are not catastrophic, but typically indicate that the 
solver is having difficulties due to numerical issues or nonmonotonicity in 
the model.

\begin{table}
\centering
\caption{Linear Solver Codes}\label{tab:7.5}
\begin{tabular}{lp{12.0cm}}
Code & Meaning\\
\hline
C & A cycle was detected.\\
E & An error occurred in the linear solve.\\
I & The minor iteration limit was reached.\\
N & The basis became singular. \\
R & An unbounded ray was encountered. \\
S & The linear subproblem was solved. \\
T & Failed to remain within tolerance after factorization was performed. \\
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Step Type Codes}\label{tab:8}
\begin{tabular}{lp{12.0cm}}
Code & Meaning\\
\hline
B & A Backtracking search was performed from the current iterate to the Newton 
    point in order to obtain sufficient decrease in the merit function. \\
D & The step was accepted because the Distance between the current iterate and
    the Newton point was small. \\
G & A Gradient step was performed. \\
I & Initial information concerning the problem is displayed. \\
M & The step was accepted because the Merit function value is smaller than
    the nonmonotone reference value. \\
O & A step that satisfies both the distance and merit function tests. \\
R & A Restart was carried out. \\
W & A Watchdog step was performed in which we returned to the last point 
    encountered with a better merit function value than the nonmonotone
    reference value (M, O, or B step), regenerated the Newton point, and 
    performed a backtracking search. \\
\end{tabular}
\end{table}

\subsubsection{Minor Iteration Log}
If more than 500 pivots are performed in the linear solver, a minor log is 
output that gives more details of the status of these pivots.  A listing 
from \texttt{transmcp} model follows, where we have set the
\texttt{output\_minor\_iteration\_frequency} option to 1.
\single\begin{verbatim}
     Minor Iteration Log
     minor      t          z     w     v   art ckpts  enter    leave
         1  4.2538e-01     8     2     0     0     0 t[    0] z[   11]
         2  9.0823e-01     8     2     0     0     0 w[   11] w[   10]
         3  1.0000e+00     9     2     0     0     0 z[   10] t[    0]
\end{verbatim}
\double
Note that \texttt{t} is a parameter that goes from zero to 1 for normal 
starts in the pivotal code.  When the parameter reaches 1, we are at a 
solution to the subproblem.  The \texttt{t} column gives the current value for
this parameter.  The next columns report the current number of problem 
variables $z$ and slacks corresponding to variables at lower 
bound $w$ and at upper bound $v$.  Artificial variables are also noted in the 
minor log.  Checkpoints are times where the basis matrix is refactored.
The number of checkpoints is indicated in the \texttt{ckpts} column.  Finally, 
the minor iteration log displays the entering and leaving variables during
the pivot sequence.

\subsubsection{Restart Log}\label{sec:restarts}
The PATH 4.x code attempts to fully utilize the resources provided by the 
modeler to solve the problem.  The algorithm tries to determine when a 
stationary point of the residual function has been encountered.  
When such a point is encountered and no progress is being made, the code
is restarted from the initial point supplied with a different set of options.  
These restarts give the flexibility to change the algorithm behavior in the 
hopes that the modified algorithm leads to a solution. The ordering and 
nature of the restarts were determined by empirical evidence based upon tests 
performed on real-world problems.  

The exact options set during the restart are given in the restart
log, part of which is reproduced below.
\single\begin{verbatim}
     Restart Log
     proximal_perturbation 0
     crash_method none
     crash_perturb yes
     nms_initial_reference_factor 2
     proximal_perturbation 1.0000e-01
\end{verbatim}
\double
If a particular problem solves under a restart, a modeler can
circumvent the wasted computation by setting the appropriate options
as shown in the log.  Note that sometimes an option is repeated in
this log.  In this case, the last value set is used.

\subsubsection{Solution Log}
A solution report is now given by the algorithm for the point returned.
The first component, under the heading ``FINAL STATISTICS'',
evaluates several different merit functions.  Next, a display of some 
statistics concerning the final point is given.  This report can be used 
detect problems with the model and solution as detailed in 
\chapref{diagnostic}.

At the end of the log file, summary information regarding the
algorithm's performance is given.  The string ``** EXIT - solution
found'' is an indication that PATH solved the problem.  Any other
EXIT string indicates termination at a point that may not be a
solution.

\subsection{Options Summary}\label{sec:path-options}
The behavior of PATH 4.x can be modified by setting particular
options.  We give a list of the available options along with their default 
and meaning in \tabref{6} and \tabref{6a}.  Note that only the 
first three characters of every word are significant.

\begin{table}
\centering\small
\caption{PATH 4.x Options}\label{tab:6}
\begin{tabular}{llp{7.25cm}}
Option & Default & Explanation \\
\hline
convergence\_tolerance & 1e-6 & Stopping criterion\\
crash\_iteration\_limit & 50 & Maximum iterations allowed in crash\\
crash\_merit\_function & fischer & Merit function used in crash method\\
crash\_method & pnewton & pnewton or none\\
crash\_minimum\_dimension & 1 & Minimum problem dimension needed 
to use crash technique\\
crash\_nbchange\_limit & 1 & Number of changes to the basis allowed\\
crash\_perturb & yes & Perturb the problem using pnewton crash\\
crash\_searchtype & line & Searchtype to use in the crash method\\
cumulative\_iteration\_limit & 10000 & Maximum minor iterations allowed\\
gradient\_searchtype & arc &  Searchtype to use when a gradient step is taken\\
gradient\_step\_limit & 5 & Maximum number of gradient step allowed before
restarting\\
interrupt\_limit & 5 & CTRL-C's required before killing job\\
lemke\_start & automatic & Frequency of lemke starts (automatic, 
first, always) \\
major\_iteration\_limit & 500 & Maximum major iterations allowed\\
merit\_function & fischer & Merit function to use (normal or fischer) \\
minor\_iteration\_limit & 1000 & Maximum minor iterations allowed in each 
major iteration\\
nms & yes & Allow line searching, watchdoging, and nonmonotone descent\\
nms\_initial\_reference\_factor & 20 & Controls size of initial reference 
value \\
nms\_maximum\_watchdogs & 5 & Maximum number of watchdog steps allowed\\
nms\_memory\_size & 10 & Number of reference values kept\\
nms\_mstep\_frequency & 10 & Frequency at which m steps are performed\\
nms\_searchtype & line & Search type to use (line, or arc)\\
output & yes & Turn output on or off.  If output is turned off, selected
parts can be turned back on. \\
output\_crash\_iterations & yes & Output information on crash 
iterations\\
output\_crash\_iterations\_frequency & 1 & Frequency at which crash 
iteration log is printed\\
output\_errors & yes & Output error messages\\
output\_factorization\_singularities & yes & Output linearly dependent columns
determined by factorization\\
output\_final\_degeneracy\_statistics & no & Print information regarding degeneracy at the solution\\
output\_final\_point & no & Output final point returned from PATH\\
output\_final\_point\_statistics & yes & Output information about the point, function, and Jacobian at the final point\\
\end{tabular}
\end{table}

\begin{table}
\centering\small
\caption{PATH 4.x Options (cont.)}\label{tab:6a}
\begin{tabular}{llp{7.25cm}}
Option & Default & Explanation \\
\hline
output\_final\_scaling\_statistics & no & Display matrix norms on the Jacobian at the final point\\
output\_final\_statistics & yes & Output evaluation of available merit 
functions at the final point\\
output\_final\_summary & yes & Output summary information\\
output\_initial\_point & no & Output initial point given to PATH\\
output\_initial\_point\_statistics & yes & Output information about the point, function, and Jacobian at the initial point\\
output\_initial\_scaling\_statistics & yes & Display matrix norms on the Jacobian at the initial point\\
output\_initial\_statistics & no & Output evaluation of available merit 
functions at the initial point \\
output\_linear\_model & no &  Output linear model each major iteration\\
output\_major\_iterations & yes & Output information on major iterations\\
output\_major\_iterations\_frequency & 1 & Frequency at which major 
iteration log is printed\\
output\_minor\_iterations & yes & Output information on minor iterations\\
output\_minor\_iterations\_frequency & 500 & Frequency at which minor 
iteration log is printed\\
output\_model\_statistics & yes & Turns on or off printing of all the 
statistics generated about the model\\
output\_options & no & Output all options and their values\\
output\_preprocess & yes & Output preprocessing information\\
output\_restart\_log & yes & Output options during restarts\\
output\_warnings & no & Output warning messages\\
preprocess & yes & Attempt to preprocess the model\\
proximal\_perturbation & 0 & Initial perturbation\\
restart\_limit & 3 & Maximum number of restarts (0 - 3)\\
return\_best\_point & yes & Return the best point encountered or the absolute 
last iterate\\
time\_limit & 3600 & Maximum number of seconds algorithm is allowed to run\\
\end{tabular}
\end{table}

The {\tt proximal\_perturbation} option can be used to overcome
problems with a singular basis matrix.  The results in a proximal
point perturbation \cite{billups:algorithms,billups.ferris:qpcomp} being
applied.  Linearly dependent columns can be output with the 
{\tt output\_factorization\_singularities} option.  For more information on 
the problems caused by singularity, we refer the reader 
to \chapref{diagnostic}.

As a special case, PATH can emulate Lemke's method 
\cite{cottle.dantzig:complementary,lemke.howson:equilibrium}
for LCP with the following options:
\single\begin{verbatim}
     crash_method none;
     crash_perturb no;
     major_iteration_limit 1;
     lemke_start first;
     nms no;
\end{verbatim}
\double
If instead, PATH is to imitate the Successive Linear Complementarity
method (SLCP, often called the Josephy-Newton method) 
\cite{josephy:newtons*1,mathiesen:algorithm,mathiesen:computation} 
for MCP with an Armijo style linesearch on the normal map residual, then 
the options to use are:
\single\begin{verbatim}
     crash_method none;
     crash_perturb no;
     lemke_start always;
     nms_initial_reference_factor 1;
     nms_memory size 1;
     nms_mstep_frequency 1;
     nms_searchtype line;
     merit_function normal;
\end{verbatim}
\double
Note that the combination of {\tt nms\_memory\_size 1} and 
{\tt nms\_initial\_reference\_factor 1} turns the nonmonotone 
linesearch off.  The {\tt nms\_mstep\_frequency 1} option disables
watchdoging \cite{chamberlain.powell.ea:watchdog}, while
{\tt nms\_searchtype line} forces PATH to search the line segment 
between the initial point and the solution to the linear model.
Finally, {\tt merit\_function normal} tells PATH to use the normal 
map for calculating the residual.

\section{Summary}
This chapter developed theory for the globalization of successive 
linearization algorithms using the Fischer-Burmeister merit function.
Algorithm GDM was shown to be well-defined for arbitrary complementarity 
problems, as well as globally and locally fast convergent under a strong 
regularity assumption.  This algorithm forms the basis for PATH 4.x which
was shown to be more robust than previous versions of the code.  In 
particular, PATH 4.x has been able to solve previously unsolvable models, 
such as the \texttt{romer} and \texttt{endog} models in MCPLIB. Finally, we 
documented the output generated by the PATH 4.x code and the options available
to users.

We remark that PATH 2.9 was used for comparative purposes in this chapter as 
it was the last release of PATH written by Steven Dirkse and Michael Ferris.
This code is based on the nonsmooth Newton method in 
\cite{ralph:global,robinson:newtons} and is documented in
\cite{dirkse:robust,dirkse.ferris:path,dirkse.ferris:pathsearch,dirkse.ferris:crash}.

The PATH 3.x series used PATH 2.9 as a template, but was completely 
rewritten.  The main modification mades were the addition of artificial 
variables, the cycle detection rules, the EXPAND pivot rules, and an 
interpolating pathsearch.  Diagnostic information and restarts were also 
introduced in the PATH 3.x series.

The PATH 4.x series is a major revision of PATH 3.x that implements 
Algorithm GDM and conforms to the solver interface in \chapref{interfaces}.  
We modified all of the searches (arc, line, and path) to 
use the Fischer-Burmeister function.  Additional enhancements made in PATH 4.x 
are the correction step used at a the solution to the linear model, updates to 
the nonmonotone search for gradient steps, and preprocessing.  We further
note that the linear solver from PATH 4.x was reused in a predictor-corrector 
method \cite{ferris.munson.ea:homotopy}.

We finally remark that an earlier version of \secref{theory} appeared in 
\cite{ferris.kanzow.ea:feasible} and parts of \secref{implementation} and 
\secref{user} appeared in 
\cite{ferris.munson:case,ferris.munson:interfaces,ferris.munson:complementarity}.

% Current research focusses on the linear solver (adaptive pivot rules, 
% one-pass pivot rules, resets, improved cycle detection and corrector step) 
% and nonlinear code (aggregations in preprocessor, interval evaluations).
