\chapter{PCPATH}
\section{Introduction}
\label{sec:intro}

The complementarity problem arises in many different applications 
\cite{ferris.pang:engineering}.
The original source of these problems were the optimality conditions
of linear, quadratic, and nonlinear programs.  Since that time, many
other applications, including game theory, economic equilibria, and
structure design/failure have been postulated and solved as
complementarity problems.   In this paper, we will
focus on the mixed (nonlinear) complementarity problem of finding 
$z \in [l,u]$ such that for all $i=1,\ldots,n$,
\[
\begin{array}{rccl}
f_i(z) &\ge  & 0 & \mbox{if } z_i = l_i % < u_i,
\\
f_i(z) & \le & 0 & \mbox{if } z_i = u_i % > l_i,
\\
f_i(z) & = & 0 & \mbox{if } l_i < z_i < u_i,
\end{array}
\]
where $f:\real^n \to \real^n$ is a continuously differentiable 
function, $l_i$ and $u_i$ are %possibly infinite valued scalars satisfying
fixed, possibly infinite numbers satisfying
\[
-\infty \;\le\; l_i \;<\; u_i \;\le\; +\infty
\]
and $[l,u]$ denotes the set of $z \in \real^n$ such that
$l_i\le z_i \le u_i$ for all $i$.  We call this problem MCP$(f,l,u)$ or 
MCP for short.

Many monotone equilibrium problems and convex optimization problems
whose constraints $C$ are {\em not} simple boxes $[l,u]$ can 
also be reformulated 
as MCPs.  For example,  consider the nonlinear program
\[
\min \phi(z) \;\mbox{ subject to }\; Az = b, \; g_j(z) \le 0 \mbox{ for }
j=1, \ldots, q
\]
where $f$ and each $g_j$ are real, differentiable,
convex functions on $\real^n$, and $A \in \real^{p\times n}$, $b \in
\real^p$.
If the Slater constraint qualification holds,
namely there exists $\hatz$ such that $A\hatz = b$ and $g_j(\hatz) < 0$
for each $j$, then finding a minimizer $z$, global or local, of this
nonlinear program is equivalent to solving its stationarity
or Karush-Kuhn-Tucker conditions; see \cite{mangasarian:nonlinear}.
The KKT conditions form an MCP in $z$ and auxiliary variables 
$\mu \in \real^p$ and $\lambda \in \real^q$
called (KKT or Lagrange) multipliers.
Define an MCP function
$f:\real^{n+p+q} \to \real^{n+p+q}$ by %of $(z,\mu,\lambda)$ by
\begin{eqnarray*}
f(x,\mu,\lambda) &=& (\nabla \phi(x) \,-\, A^T \eta \,-\, 
        \nabla g(x)^T\lambda,
      %\sum_j \nabla g_j(x)\lambda_j, 
\, Ax-b, \, g(x))
%\\
%l & = & (-\infty^n, \, -\infty^p, \, 0^q)
%\\
%u & = & (infty^n, \, \infty^p, \, \infty^q)
\end{eqnarray*}
where $\nabla \phi(x)$ is the gradient vector of $\phi$ at $x$,
$g(x) = (g_1(x), \ldots, g_q(x))$ and $\nabla g(x)$
is the $q\times n$ Jacobian matrix of $g$ at $x$.
Also let the box $[l,u]$ be $\real^n\times \real^p \times \real^q_-$.
Then the associated MCP is precisely the KKT conditions of the
nonlinear program, hence is equivalent to the nonlinear program.

Our method will be based on a reformulation of the mixed
complementarity problem.  To explain this reformulation, let us first
consider a nonempty polyhedral subset of $\real^n$, $C$.
Associated with each of the faces of $C$ is a (full-dimensional)
polyhedral set, and the collection of these polyhedra comprise a
piecewise linear manifold of $\real^n$ called the {\em normal
manifold}, denoted $\calNC$. These polyhedra are called the {\em
cells} of $\calNC$; a full description along with important
properties is given in \cite{robinson:normal} (see also
\cite{ralph:branching} for further investigation).
For example, when $C = \real^n_+$, that is $C =[l,u]$ with each $l_i=0$
and $u_i=\infty$, the cells of $\calNC$ are the orthants
of $\real^n$.  

We denote by $\pic{\cdot}$ the Euclidean projection mapping onto the 
set $C$; thus for the box $C=[l,u]$, $\pic{x}$ is the vector
whose $i$th component is $l_i$ if $x_i\le l_i$,
$x_i$ if $l_i\le x_i \le u_i$, and $u_i$ otherwise.
The normal map 
\cite{robinson:normal} induced by $(f,C)$ is the function 
$\FC:\real^n \to \real^n$ given by
\[
\FC(x) \;=\; f(\pic{x}) + x - \pic{x}.
\]
It is well known that the projection map $\pic{x}$
and hence the normal map $\FC(x) = f(\pic{x}) + x - \pic{x}$
are smooth in the interior of each of the cells of $\calNC$.
The key point is that for $C = [l,u]$, if $x$ satisfies 
$\FC(x) = 0$, then $z = \pic{x}$ solves MCP.  
Furthermore, if $z$ solves MCP then $x = z - f(z)$ is a zero of $\FC$.

We propose and implement a homotopy (or continuation) method for
finding zeroes of the normal map $\FC(x)$ when $C = [\ell,u]$
based on classical predictor-corrector
ideas. Previously proposed algorithms
\cite{sellami.robinson:homotopies,sellami:continuation,sellami.robinson:implementation} 
use linear
approximations to single (smooth) pieces of a piecewise smooth
mapping, see also \cite{kojima:recent,kojima.hirabayashi:homotopy}.  
The method proposed in this 
paper is based on piecewise
linear approximations of a piecewise smooth mapping.  Both our method
and that of \cite{sellami.robinson:implementation} are based upon the
theoretical foundations developed in
\cite{alexander:topological,alexander.kellog.ea:piecewise,alexander.li.ea:piecewise}.

%%% need to mention Billups work.  I made up the following.
%%% Can you check this against his paper?
An alternative homotopy approach investigated recently \cite{billups.speight.ea:nonmonotone}
is to form a homotopy that is smooth except in the limit as the
homotopy parameter reaches its final value.  This makes application
of smooth homotopy codes attractive
though some care near the end of the homotopy path, when the problem
is tending toward nonsmoothness, may be needed.

\secref{background} gives some background material on homotopy
methods.  \secref{homotopy} outlines the theoretical underpinnings of
our method.  \secref{implementation} describes the implementation of a
large scale code based on this theory and \secref{results2} gives some
computational results on problems from MCPLIB
\cite{dirkse.ferris:mcplib}.

\section{Background}
\label{sec:background}

We give the notation and background for homotopy methods,
in particular for normal maps induced by smooth functions
and polyhedral convex sets.

\subsection{Homotopy mappings and paths}
Continuation or homotopy methods are a technique to trace the zeroes of a 
homotopy mapping, starting from an easily constructed zero and moving towards
the solution to the problem of interest.  To describe such methods, we
first need to understand the homotopy mapping.

For a fixed vector $a \in \real^n$,
we construct a homotopy mapping $\HC^a:\real^{n+1}\to \real$ 
from the function $x\mapsto x-a$ to the function $\FC$ by interpolation:
for $(x,t) \in \real^{n}\times [0,1]$,
\[
\HC^a(x,t) \;=\; (1-t)(x-a) + t \FC(x).
\]
Note that $x=a$ is the unique solution of $\HC^a(x,0) = 0$,
and $\HC^a(x, 1) = 0$ if and only if $x$ is a zero of the
normal map.  Hence if we can calculate the endpoint $x^a(1)$
of the homotopy path, then we have solved the original
problem.  
The idea of a continuation method is to analytically or numerically
determine a path of solutions $x^a(t)$ of the equation
$\HC^a(\cdot,t) = 0$ from $t=0$, where $x^a(0) = a$ is the unique solution,
to $t=1$, where as mentioned, $x^a(1)$ is a zero of the normal map.

A standard list of properties of the homotopy path $x^a(t)$
follows.
We only state these for the case of $\HC^a$ defined
above, though the class of homotopies with these properties
is considerably larger
\cite{alexander.li.ea:piecewise} as we explain after
Theorem~\ref{thm:path existence}.
Parts~1 and~2 below summarize results that are due, in essence,
to \cite{alexander.li.ea:piecewise}, while part~3 is due to
\cite{sellami:continuation}; see the proof for details.
Parts~4 and~5 seem to be new.  Part~5 refers to subanalytic functions
\cite{lojasiewicz:ensembles,hironaka:introduction}
which we define,
for completeness, following the statement of the theorem.
However we
refer to \cite[Section 5.2]{facchinei.pang:total}
for a succinct introduction to the properties of subanalytic
sets and functions, and additional references.
Also, the notation $t \to 1^-$ means $t \to 1$, $t < 1$.

\begin{theorem}\label{thm:path existence}
Let $f:\real^n \to \real^n$ be {\rm{C$^1$}} and $C$ be a nonempty polyhedral
convex set in $\real^n$.  For almost all $a\in \real^n$:
\begin{enumerate}
\item \label{thm:path:1}
There exist $T > 0$  and a piecewise smooth
path $x^a:[0,T) \to \real^n$ such that $x^a(0) = a$,
and $\{(x^a(t),t) \,;\, t \in [0,T)\}$
is a connected component of the set 
$\{(x,t)\,:\, t \in [0,T),\, \HC^a(x,t) = 0 \}$.

\item \label{thm:path:2} The function $x^a(\cdot)$ is nondegenerate in
the sense that for each $t \in [0,T)$, first, $x^a(t)$ lies either in the 
interior of a unique cell of the normal manifold $\calNC$ or in the
relative interior of a facet (an $(n-1)$-dimensional face of a cell)
that is the intersection
of two cells of $\calNC$; and, second,
if $x^a(t)$ lies in a cell $\sigma$ and 
$g_\sigma$ denotes the smooth mapping that represents $\HC^a$ on $\sigma$,
then the Jacobian (derivative) matrix  $\nabla g_\sigma(x^a(t),t)$
has full rank.

\item \label{thm:path:4} If $f$ is {\rm C}$^2$ and 
there exists a point $u \in C$ and a scalar $\rho > \norm{u}$
such that for each $c \in C$ with $\norm{c} = \rho$ we have
$\ip{f(c)}{c-u} \ge 0$, then 
$\{x^a(t):t\in [0,1)\}$ is bounded and its closure 
contains a point $\hatx$ with $\FC(\hatx) = 0$.

\item\label{thm:path:4a}
If $f$ is {\rm C}$^2$, $T \ge 1$, and there exists a limit point $\hatx$
of $p(t)$ as $t\to 1^-$ such that the directional derivative
$\FC'(\hatx;\cdot)$ is invertible, then the arc length
of the path $x^a:[0,1) \to \real^n$ is finite and 
$x^a(t) \to \hatx$ as $t \to 1^-$.

\item \label{thm:path:3} If $T \ge 1$ and $\hatx$ is a limit point 
of $x^a(t)$ as $t \to 1^-$,
then $\FC(\hatx)= 0$.
If, in addition, either $T > 1$, or $f$ is subanalytic
then $x^a(t) \to \hatx$ as $t \to 1^-$.
\end{enumerate}
\end{theorem}

\begin{pf}
Parts~1 and~2 are essentially due to \cite[Theorem 1]{alexander.li.ea:piecewise}.
The fact that the graph of $x^a$ is a connected
component of $\{(x,t)\,:\, t \in [0,T),\, \HC^a(x,t) = 0 \}$
is not stated but can be easily deduced from
the discussion on transversality prior to \cite[Theorem 1]{alexander.li.ea:piecewise}.
See also \cite[Theorem~1]{sellami.robinson:implementation}
and subsequent discussion of these transversality properties.
Part~\ref{thm:path:4} is quoted from \cite[Proposition~5.1]{sellami:homotopy}.

For part~\ref{thm:path:4a}, suppose $\hatx$ is a limit point
of $x^a(t)$ as $t\to 1^-$ such that $\FC'(\hatx;\cdot)$ is
invertible.  The implicit function theorem of \cite{robinson:implicit-function} 
says that the path $x^a(t)$ is locally uniquely defined for $t$ near $1^-$:
for some $\epsilon > 0$ and all $t$ near 1 with $t< 1$, the homotopy equation 
$\HC^a(\cdot,t) = 0$ has a unique solution $x$ within distance $\epsilon$ of $\hatx$.
It follows that the path $\{(x^a(t),t) \,:\, t \in [0,1)\}$
is bounded and the only limit point of $x^a(t)$ as $t \to 1^-$ is
$\hatx$ itself; i.e.~$x^a(t)$ converges to $\hatx$ as $t \to 1^-$.
Now invertibility of $\FC'(\hatx;\cdot)$ is 
equivalent to other properties such as strong regularity of the
associated variational inequality at
$\pic{\hatx}$, \cite[Definition~2.2]{sellami.robinson:implementation},
and coherent orientation of the normal map $\nabla f(\hatx)_{K(\hatx)}$,
\cite{robinson:normal}, where $K(\hatx)$ is the critical cone to $C$
at $\hatx$, and this normal map actually coincides with
$\FC'(\hatx;\cdot)$; see \cite[Section~5]{sellami:continuation} and
\cite[Section~2]{sellami.robinson:implementation} for discussion and
related results.
Therefore the analysis of \cite[Proposition~5.6]{sellami:continuation}
can be adapted to show finite arc-length, using boundedness
of the path and the fact that $\hatx$ is the only limit point of
$x^a(t)$ as $t\to 1^-$.

For part~\ref{thm:path:3}, it is clear by continuity of $\FC$ that $\FC(\hatx) = 0$
if $(\hatx,1)$ lies in the closure of $\{(x^a(t),t) \,:\, t \in
[0,1)\}$.  If $T > 1$, continuity of $x^a(\cdot)$ immediately yields
convergence of $x^a(t)$ to $\hatx$ as $t \to 1$.
Suppose instead that $f$ is subanalytic. 
We will show convergence of $x^a(t)$ to $\hatx$ as
$t \to 1^-$ by refining an argument in the proof
of \cite[Theorem~5.10]{facchinei.pang:total}.
All polyhedral convex sets are
subanalytic so that the projection $\pi_C$ mapping is also subanalytic
\cite[Lemma~5.9]{facchinei.pang:total}.  It follows that $\HC^a$ is also
subanalytic since the composition of subanalytic mappings is also
subanalytic, hence that $(\HC^a)^{-1}(0)$ is a subanalytic set.  Let
$P^a = \{(x^a(t),t) \,:\, t \in [0,1)\}$; from part~\ref{thm:path:1},
taking $T = 1$,
%for each $t \in [0,T)$ there is $\epsilon_t > 0$ such that 
there is a neighborhood $U$ of $P^a$ such that $U \cap (\HC^a)^{-1}(0)
= P^a$.  Without loss of generality assume $U$ is subanalytic; so then
$P^a$, the intersection of subanalytic sets, is also subanalytic.
This shows that $x^a:[0,T) \to \real^n$ is subanalytic.
The lemma in the Appendix concludes the proof.
\end{pf}

With regard to part~\ref{thm:path:3}: a set $S$ in $\real^n$ is {\em
subanalytic} if there exists a semianalytic set $T$ in a higher
dimensional space $\real^{N}$ such that $S = \{x \in \real^n \,:\,
(x,y) \in T\}$.  A set $T$ in $\real^N$ is {\em semianalytic} if for
each $\barz \in T$ there is a neighborhood $V$ of $\bar{z}$ such that
$T \cap V$ can be written as the finite union of sets of the form
%
\[ V \cap \{z \in \real^N \,:\, f_i(z) = 0, i =1,\ldots, I; \, g_j(z)
< 0, j = 1,\ldots, J \} \]
%
where $I$ and $J$ are nonnegative integers, and each $f_i$ and $g_j$
is a real analytic mapping on $\real^N$.  A function is subanalytic if
its graph is a subanalytic set.  Thus the class of subanalytic
functions is rather broad.

Much of Theorem~\ref{thm:path existence} holds
if we replace the homotopy mapping $\HC^a$ by one of the form 
\[ 
(1-t)\phi(a,x) + t \FC(x)
\]
where $\phi:\real^{m+n}\to\real^n$ is a {\em sufficient}
mapping \cite{alexander.li.ea:piecewise}, that is
$\phi$ is smooth 
such that its partial Jacobian matrix with respect
to $a \in \real^m$, $\nabla_a\phi(a,x)$, has full rank for all $(a,x)
\in \real^m\times\real^n$. (For example, the function $\phi(a,x) =x-a$ 
used in the definition of $\HC^a$ is sufficient.)  To be precise, 
in order to generalize parts~\ref{thm:path:1}, \ref{thm:path:2} 
and~\ref{thm:path:3} of Theorem~\ref{thm:path existence} using
a homotopy with a sufficient mapping as described,
we only need to replace the qualifier
$a \in \real^n$ by $a \in \real^m$ and, in the statement of
part~\ref{thm:path:3}, 
% subanalyticity, to require subanalyticity of $\phi$.
the condition ``$f$ is subanalytic'' by ``$f$ and $\phi$ are
subanalytic''.



For future reference, we define the homotopy path as the set
\[
P^a \;=\; \{(x^a(t),t))
\,:\, t \in [0,T)\}
\]
where $T$ is the maximum value in $(0,1]$ such that part~1 of
Theorem~\ref{thm:path existence} holds.

\subsection{Predictor-corrector methods}
Predictor-corrector methods are a class of 
algorithms that attempt to numerically traverse
the homotopy path by a sequence of 
predictor-corrector iterations. Suppose for the moment that
$\HC^a$ is smooth, e.g.~$C = \real^n$ and $\FC = f$.
At step $k$ we are given an iterate $(\xk,\tk)$ that is approximately
on the path, that is $\xk \approx x^a(\tk)$.
The predictor step identifies a nonzero vector 
$d = (d_x,d_t) \in \real^{n}\times \real^1$ 
tangent to the path at the current point % $(x^a(t),t)$,
and estimates a new point further along the path
$(\xkd,\tkd) = (\xk,\tk) +hd$ for some small $h >0$.
The direction $d$ is found as a vector in the kernel
of the Jacobian matrix $\nabla \HC^a(\xk,\tk)$ that
maintains the correct orientation, that is we want
to move tangent to the curve in the direction that
increases arc length.
The corrector step then
tries to identify a point $(x^a(s),s)$ on
the path near to $(\xkd,\tkd)$; this is used to define
the next iterate, $(\xkp,\tkp) = (x^a(s),s)$.
The corrector step is usually carried out by a version
of Newton's method that uses the Moore-Penrose inverse
of the $n \times (n+1)$ matrix $\nabla \HC^a(x,t)$
starting with $(x,t) = (\xkd,\tkd)$ and proceeding
until $\HC^a(x,t)$ is approximately zero;
see \cite{allgower.georg:numerical}.

When $\HC^a$ is piecewise smooth, as it is when $f$
is C$^1$ and $C$ is polyhedral convex,
it is possible to mimic the predictor and corrector steps taken
in the smooth case by either staying within one cell, or identifying
an adjacent cell and moving into it.  This means that
at each predictor step a cell $\sigma_k$ containing the current
iterate $(\xk,\tk)$ is identified, and the predictor direction $d$
is chosen using the Jacobian matrix of the smooth mapping
corresponding to $\HC^a$ on $\sigma^k$.
The Moore-Penrose-Newton iteration is undertaken in a similar way.
This is the line of development in \cite{sellami.robinson:implementation}.
We note that these algorithms, and our proposed method to follow,
include some nontrivial technical details that are necessary to make it 
theoretically viable. In particular,
it is assumed that the chosen starting point $a$ yields a path 
$x^a(t)$ that is nondegenerate as described
in Theorem~\ref{thm:path existence}, part \ref{thm:path:2}.
The implementation must take special care in the near-degenerate case 
for which the path
passes through or near a point that is contained in three or more cells.

We prefer to use a piecewise linear (PL) predictor approximation,
that is to generate a piecewise linear tangent path to
the piecewise smooth homotopy path.  An advantage of the PL predictor is that,
in principle, it allows the near degenerate situation to be handled 
in the same way as the nondegenerate case. 
Another advantage of the PL tangent path approach is that we
can use the code developed over several versions of the PATH 
solver \cite{dirkse.ferris:path} for MCP; we rely on this code to handle 
numerical degeneracy in the PL path for example.

Consistency would dictate that each step of the corrector also be 
taken with respect to a PL model of the actual path. However for 
simplicity of implementation we
have chosen a method more like that of \cite{sellami.robinson:implementation}:
we pick a cell containing the current point and apply a Moore-Penrose-Newton 
iteration as a heuristic for decreasing the distance to the path.
Details of both the predictor and corrector step are given below.

\section{Homotopy method}
\label{sec:homotopy}


\subsection{Predictor step}
\label{sec:predictor}
To describe the predictor step, 
we need to define a piecewise linear path that
is tangent to the piecewise smooth homotopy path $P^a$ at
a given point.  To do this, we need to describe a linearization of $\HC^a$.

For convenience in approximating $\HC^a$, we write it in another way.
Let $H^a:\real^{n+1} \to \real^n$ be the C$^1$ homotopy from the 
mapping $z \mapsto z-a$ to $f(z)$,
\[
{H^a}(z,t) \;=\; (1-t)(z-a) + t f(z),
\]
and observe that
\[
\HC^a(x,t) \;=\; H^a(\pic{x},t) \,+\, x-\pic{x}.
\]

We define the linearization of the smooth mapping $H^a$ at any 
point $(z,t) \in \real^{n+1}$
using the first two terms of the Taylor's series expansion:
for $(z',t') \in \real^{n+1}$, 
\[
\LH^a(z,t)(z',t') \;=\; H^a(z,t) 
                     \,+\, 
%\nabla H^a(z,t)(z'-z,t'-t),
%[(1-t)I + t \nabla f(z) \;\; f(z)-(z-a)].
[(1-t)I + t \nabla f(z)](z'-z) \,+\, [f(z)-(z-a)](t'-t)
.
\]
%where $\nabla H^a(z,t) \in \real^{n\times(n+1)}$ is the Jacobian of
%$H^a$ at $(z,t)$.  
%
%We have
%\[
%\nabla H^a(z,t) \;=\; [\nabla_zH^a(z,t) \;\; \nabla_tH^a(z,t)]
%\]
%where $\nabla_z,\nabla_t$ are the partial derivatives with respect
%to $z,t$ respectively.  Thus
%\[\begin{array}{rcccl} %\beqas
%\nabla_zH^a(z,t) & = &
%(1-t)I + t \nabla f(z) & \in & \real^{n\times n} \\
%          \nabla_tH^a(z,t) & = & f(z) -(z-a) & \in & \real^n
%\end{array}\] %\eeqas
%where $I$ is the $n\times n$ identity matrix.

Now we define the linearization of the piecewise smooth mapping 
$\HC^a$ at $(x,t)$ as a function of 
$(x',t') \in \real^{n+1}$,
\[
\calL \HC^a(x,t)(x',t') \;=\;  
\calL H^a(\pic{x},t)(\pic{x^\prime},t') \,+\, x' - \pic{x'}
.
\]
This is a piecewise linear mapping in $(x',t')$.
It is a ``point based approximation'' to $\HC^a$ at $(x,t)$ in the 
terminology of \cite{robinson:newtons}.

By the ``piecewise linear path that
is tangent to the piecewise smooth homotopy path $P^a$ at
$(\xk,\tk)$'' we mean the set
\begin{equation} \label{eqn:predictor}
P^a_k \;=\; \{(x,t) \,:\, t \mbox{ near } \tk,\, \LHC^a(\xk,\tk)(x,t) = 0 \}.
\end{equation}
Note $(\xk,\tk) \in P^a_k$.
Given a step size $h_k > 0$ our task is to partially 
construct $P^a_k$, starting from $(\xk,\tk)$ 
and moving
in the direction that is associated with increasing
the arc-length of the actual homotopy path $P^a$. This ``one-sided''
path is denoted $Q^a_k$.
We keep moving along $Q^a_k$ until we 
determine a point $(x',t')$ at distance $h_k$ from $(\xk,\tk)$.
The new point $(x',t')$
is our prediction for the next point on the path $P^a$.  We need
a corrector scheme to move from $(x',t')$ to a nearby
point in $P^a$; this will be discussed in \secref{corrector}.

For the subsequent discussion we assume for all $t \in [0,1)$ that
either $x^a(t)$ lies in the interior of a cell of $\calNC$
or in the relative interior of one of its facets,
and that the Jacobian matrices at $(x^a(t),t)$
of the smooth functions representing $\HC^a$ on these cells have full rank.
These properties hold for almost all $a$ by Theorem~\ref{thm:path existence}.

We note that the full rank of the Jacobians
of $\HC^a$ near $(\xk,\tk)$, together with nondegeneracy of $\xk$,
is enough to show that the set $P^a_k$
does indeed define a piecewise linear path, i.e.~a one-dimensional
piecewise affine manifold, for $t$ near $\tk$.
This PL path is tangent to the piecewise smooth homotopy
path $P^a$ at $(\xk,\tk)$ if $(\xk,\tk) \in P^a$.
It will be necessary to generate $Q^a_k \subset P^a_k$ starting from 
$(\xk,\tk)$  and moving in one of the two possible directions as determined
by the orientation calculation presented next.
In practice we allow for $(\xk,\tk)$ to be near rather than in
$P^a$ --- c.f.~the tolerance $\epsc$ used in the corrector
step to follow --- so that $P^a_k$ is approximately tangent to $P^a$.
Such details, including the invariance of orientation on the homotopy path, are
spelled out in \cite{alexander:topological,alexander.kellog.ea:piecewise,alexander.li.ea:piecewise,allgower.georg:numerical}.

\subsubsection{Orientation of the path}
\label{sec:orientation}
An orientation parameter $\eta = \pm 1$ is used to
decide which way to go on the path $P^a_k$ when starting
at $(\xk,\tk)$.
The orientation parameter is defined prior to
the first predictor step of the homotopy method as
\[\eta \;=\;
\sgn  \det\left[\begin{array}{cc} \nabla\HC^a(a,0) \\ (d^0)^T
\end{array}\right] ,
\]
where $d^0 = (-f(a),1)$.
Note $\HC^a$ is differentiable at any $(x,0)$.
Initially $t^0=0$, we have $x^0 = a = x^a(0)$, and
$\nabla\HC^a(a,0) = [I \,  f(a)]$,
so that the chosen vector $d^0$ is in the kernel of $\nabla\HC^a(x^0,t^0)$
as in the smooth case.  
Therefore
\begin{eqnarray*}
\eta &=& \sgn \det \left[\begin{array}{cc} I & f(a) \\ -f(a)^T & 1
\end{array}\right]
\\ & = & \sgn\det 
\left(
\left[\begin{array}{cc} I & 0 \\ 
                        -f(a)^T & 1 + \norm{f(a)}^2
      \end{array}\right]
\left[\begin{array}{cc} I & f(a) \\ 0 & 1\end{array}\right]\right)
\;=\; 1 ,
\end{eqnarray*}
where we are using the fact that the determinant of a product of two
matrices is the product of the determinants of the matrices.
Writing $d^0 = (d_x,d_t) \in \real^n\times \real$,
we have chosen $d_t =1 >0$ 
since initially we plan to 
move along the path $P^a$ by increasing $t$,
i.e.~our first prediction will be a point $(a,0) + h_0d^0$,
for some $h_0>0$, which we believe will be close to
a point $(x^a(t),t)$ with $t \approx hd_t$.



Let $(\xk,\tk) \in P^a$,
$\sigma$ be a cell containing $\xk$,
and $\LHC^a(\xk,\tk)$ be represented on $\sigma$
by an affine map
whose Jacobian $J \in \real^{n\times(n+1)}$
at $(\xk,\tk)$
has full rank. So the kernel of $J$ is a one-dimensional
space.
Let $\ker J$ be spanned by a nonzero 
vector $d=(d_x,d_t) \in \real^n\times\real$,
where we also ensure $\xk + sd_x \in \sigma$ for small $s > 0$
by using $-d$ instead of $d$ if necessary
(this is possible since $\xk$ is interior either to the cell $\sigma$
or to the union of two cells).
Now calculate the sign $\eta_k$ of the determinant of the
$(n+1)\times (n+1)$ matrix consisting of the matrix $J$ augmented by
the row $d^T$ (see \secref{impl-orient}). 

If $\eta_k = \eta$ then we generate the first part of
the path $Q^a_k$  starting
at $(\xk,\tk)$ and moving in the nonzero direction $d^0 = d$
specified by the (kernel of the) Jacobian of $\LHC^a(\xk,\tk)$ in $\sigma$.
Otherwise $\eta_k = -\eta$,
and we move in the ``opposite direction'' which is
either $d^0=-d$ if $\xk \in \intr\, \sigma$ or, if $\xk$
is in the intersection of $\sigma$ and another cell $\sigma'$,
a nonzero direction $d^0$ constructed
with reference to the cell $\sigma'$ instead of $\sigma_k$.
(Invariance of orientation is demonstrated, in part,
by the result that if $\sigma^0$ denotes the cell $\sigma$ or $\sigma'$ 
associated
with $d^0$, depending on which situation occurs,
and $J^0$ is the Jacobian at $(\xk,\tk)$ of
$\LHC^a(\xk,\tk)$ restricted to $\sigma^0$, then 
the matrix consisting of $J^0$ augmented with the row $(d^0)^T$
has determinantal sign equal to $\eta$.)


\subsubsection{Generating the path $Q^a_k$}
\label{sec:generating}
The PL path $Q^a_k \subset P^a_k$ is generated, 
one line segment at a
time, in the following way.
The path generation procedure is initialized
with $(\xi^0,\tau^0) = (\xk,\tk)$,
the nonzero direction $d=d^0$
as described in \secref{orientation},
and an associated cell $\sigma^0$.
In fact, at
each step we are given a point $(\xi,\tau) \in P^a_k$, a cell
$\sigma$ of the normal manifold
$\calN_C$ that contains $\xi$,
and a nonzero direction $d = (d_x,d_t) \in \real^n\times \real$ 
such that $\xi+sd_x \in \sigma$ for
all small $s > 0$. The next line segment on the path is generated
by finding the
maximum value, possibly infinite,
of $s \ge 0$ such that $\xi + sd_x \in \sigma$; denote
this value by $s'$.  If $s'$ is finite and the path $Q^a_k$
is nondegenerate (as it must be for almost all $a$) 
then the new point is $(\xi',\tau') = (\xi,\tau) + s'd \in Q^a_k$
and there is a cell $\sigma'$
distinct from $\sigma$
such that $(\xi',\tau')$ lies in the relative interior
of a common facet of $\sigma'$ and $\sigma$.
It is now possible, as in \secref{orientation},
 to choose a unit direction $d'= (d_x',d_t')$
in the kernel of the Jacobian $J'$ of $\LHC^a(\xk,\tk)$ when
restricted to $\sigma'$, such that $\xi' + sd_x' \in \sigma'$
for small $s \ge 0$. 
The remainder of the path is generated inductively.

A more formal justification of this path generating procedure
can be given in terms of \cite{eaves:short}
where it is also shown how the degenerate case can be handled.
It follows from this paper
that $Q^a_k$ is composed
of finitely many line segments, the last of which is either
a ray or contains the starting point $(\xk,\tk)$.

The path generating procedure may be terminated in
any step that generates $(d_x,d_t)$ and $s'$ from $(\xi,\tau)$ and 
$\sigma$ as follows.
Recall $h_k$ is a positive parameter bounding the distance
of points on $Q^a_k$ to the starting point $(x^k,t^k)$.
Termination occurs at
the point $(x',t') = (\xi,\tau) + s(d_x,d_t)$ if there exists
$0 < s \le s'$ such that 
any one of the following conditions holds:
\begin{enumerate}
\item[(i)] $\norm{(x',t') - (\xi^0,\tau^0)} = h_k$
\item[(ii)] $t'=1$
\item[(iii)] $t'=0$
\item[(iv)] $(x',t') = (\xi^0,\tau^0)$.
\end{enumerate}
Stopping condition (iii) is used since there is only one solution
$x=x^0=a$ to the equation $\HC^a(x,0) =0$,
and we are not interested in revisiting $(x^0,0)$.
Stopping condition~(iv) is used to prevent a kind of cycling, that is 
traveling in an endless loop back.

\subsection{Corrector step}
\label{sec:corrector}

We are given a corrector tolerance $\epsc > 0$, 
that is we expect every iterate $(\xk,\tk)$
to satisfy
\[
\norm{\HC^a(\xk,\tk)} \;\le\; \epsc.
\]
Therefore if the point $(\xkd,\tkd)$ is the result
of a predictor step from $(\xk,\tk)$,
then the role of the corrector
step is to find the next iterate $(\xkp,\tkp)$ somewhat
near $(\xkd,\tkd)$ such that 
$\norm{\HC^a(\xkp,\tkp)} \;\le\; \epsc$.
If an iterative  corrector scheme is unable to provide
such a point within $\Jc$ iterations, where $\Jc$
is a positive integer,
then the corrector step is deemed to have failed,
and a new predicted point $(\xkd,\tkd)$ must be provided.

\subsubsection{A theoretically robust corrector scheme}
The corrector step can be extended from the Moore-Penrose Newton 
method~\cite{allgower.georg:numerical} for smooth operators, which
we paraphrase as follows.
Suppose $H:\real^n \to \real^{n+k}$ is a smooth function,
and $z^0 \in \real^{n+k}$ is given; for continuation methods
it is enough to take $k =1$.  For each $j=0, 1, 2, \ldots$, 
we assume the Jacobian matrix $\nabla H(z^j) \in \real^{(n+k)\times n}$
has full rank, and define $z^{j+1}$ as the nearest point, in
Euclidean distance, to
$z^j$ 
in $\{z \,:\, H(z^j) + \nabla H(z^j)(z-z^j) = 0\}$.
This geometrically motivated algorithm
is usually stated, equivalently, as
$z^{j+1} = z^j -\nabla H(z^j)^\dagger H(z^j)$
where $\nabla H(z^j)^\dagger$ is the Moore-Penrose inverse of
$\nabla H(z^j)$.  It converges locally at a Q-quadratic rate to
a zero $z^*$ of $H$ if $\nabla H(z^*)$ has full rank
and $\nabla H$ is Lipschitz near $z^*$.

The paper
\cite{dontchev:local}
provides an extension of this algorithm, designed for
solving generalized equations, that can be easily
adapted to finding zeroes of the nonsmooth function $\HC^a$. 
This approach has also been studied in some detail
in \cite{kummer:metric}.
We give a version of the Newton method of \cite{dontchev:local} for the case
of the corrector
step. It uses the parameters $\epsc > 0$ and $\Jc \in \natnum$.


\begin{quote}
{\bf Corrector Step A}, a function of $(x',t')$.

Let $(y^0,s^0) = (\xkd,\tkd)$ and $j=0$.

While $\norm{\HC^a(y^j,s^j)} > \epsc$ and $j < \Jc$
\begin{quote}
  Find a globally nearest point $(y^{j+1},s^{j+1})$
to $(x',t')$ in the set
\begin{equation}\label{eqn:corrector}
\begin{array}{rcl}
\{ (y,s) \,:\, \LHC^a(y^j,s^j)(y,s)  =  0 \}.
\end{array}
\end{equation}


Let 
  $j = j+1$.
\end{quote}
(end While)

Let $(\yjc,\sjc) = (y^j,s^j)$.
\end{quote}
Using Theorem~\ref{thm:path existence} it follows that
for almost every $a$, if $(y^j,s^j)$ is near enough to
the path $P^a$ then the linearized set defined
in \eqnref{corrector} is a piecewise linear path at least
in a neighborhood of $(y^j,s^j)$.  Thus determining
a globally nearest point to this path (within a closed neighborhood
of $(y^j,s^j)$) is made computationally possible by examining each of
the finitely many line segments of the path.

Suppose $f$ is a C$^2$ function.
It is known \cite{dontchev:local,kummer:metric} that if 
$\LHC^a$
satisfies certain regularity conditions, then the above Newton
method converges Q-quadratically to a zero of $\HC^a$.
Suppose further that the homotopy path $P^a$ is bounded 
and nondegenerate such that, 
for some point $(\hatx,\hatt)$ in its closure, the 
partial directional derivative of $\HC^a$ with respect to $x$,
$(\HC^a)'_x(\hatx,\hatt;\cdot)$, is invertible.
Then we claim the following statement can be established
using the Newton convergence result, Theorem~\ref{thm:path existence} 
and compactness arguments: %that
there exist
$h > 0$, $\epsc > 0$, and $\Jc \in \natnum$
such that for any $t \in [0,\hatt)$ and $x$ with 
$\norm{\HC(x,t)} \le \epsc$,
if the predictor step generates $(x',t')$ from $(x,t,h)$,
then Corrector Step~A
terminates in one step 
and the distance of $(y^1,s^1)$ to the path $P^a$ is 
a small order of the distance from $(y^0,s^0)= (x',t')$ to $P^a$.

We comment that the choice of a viable $\epsc$ seems to depend on
the ``how far the path is from degenerate points'', hence
on $a$, so it is not clear how to choose $\epsc$ appropriately
other than as a small positive number.

\subsubsection{A heuristic corrector scheme}
\label{sec:heuristic}

To simplify our implementation, that is to avoid computation of a globally 
nearest point to the generally nonconvex set \eqnref{corrector},
it is easy to extend the approach of Kojima and Shindo
\cite{kojima.shindo:extensions} for piecewise smooth systems
in the same number of variables as equations
to underdetermined piecewise smooth equations.  At the $j$th iterate,
the idea is to identify any cell $\sigma_j$ containing
the current point $(y^j,s^j) \in \real^{n+1}$,
and apply one step of the Moore-Penrose Newton method
to the smooth mapping representing $\HC^a$ on $\sigma_j$.

\begin{quote}
{\bf Corrector Step B}, a function of $(x',t')$.

Let $(y^0,s^0) = (\xkd,\tkd)$ and $j=0$.


While $\norm{\HC^a(y^j,s^j)} > \epsc$ and $j < \Jc$
\begin{quote}
  Let $(y^{j+1},s^{j+1}) \,=\, (y^j,s^j) - \nabla H_{j}(y^j,s^j)^\dagger \HC^a(y^j,s^j)
  $
\\
  where
  $\sigma_j$ is a cell of $\calN_C$ containing $(y^j,s^j)$, and
  $H_{j}$ is the C$^1$ representation of $\HC^a$ on $\sigma_j$.
\\
  Let 
$j = j+1$.
\end{quote}
(end While)

Let $(\yjc,\sjc) = (y^j,s^j)$.
\end{quote}

We stress that Corrector Step~B is only a heuristic for piecewise
smooth systems because the convergence analysis is not as strong as it is for
the classical Newton's method or its extensions in 
\cite{dontchev:local,kummer:metric}.  See \cite{chen.nashed.ea:convergence}
for a useful but limited convergence theory of this method and 
some generalizations. 

\subsection{Formal Homotopy Algorithm}
As described above, 
if a predictor step produces a point
that cannot be sufficiently corrected in $\Jc$ or fewer iterations of
the corrector step, then a revised, more conservative prediction
is made, and the correction step is attempted again.  We define this 
formally using the constants $\epst \ge 0$, $ \epsf > 0$ 
and $\lambda \in (0,1)$ at the outer level of the algorithm;
$\eta \in \{\pm 1\}$ in the predictor step;
and $\epsc > 0$ %, $\gamc > 0$ 
and $\Jc \in \natnum$ in the corrector step.

\begin{quote}
{\bf Homotopy algorithm} for solving $\HC^a(x,1) = 0$.

Let $(x^0,t^0) = (a,0) \in \p$, 
$k = 0$, $\hath \in (0,1]$, $\lambda \in (0,1)$.

While $|\tk-1| > \epst$ or $\norm{\FC(x^k)} > \epsf$
\begin{quote}

Let $\hk=\hath$.

Repeat:
\begin{quote}

$(x',t',\hk) \;\longleftarrow$ {\bf Predictor step}$(\xk,\tk,\hk)$.

$({x^\dagger},{t^\dagger}) \;\longleftarrow$ {\bf Corrector step B}$(x',t')$. 

Let $h_k = \lambda h_k.$
\end{quote}
Until  $\norm{\HC(\yjc,\sjc)} \le  \epsc$.

Let $(\xkp,\tkp) = (\yjc,\sjc)$, $k = k+1$.

\end{quote}
(end While)
\end{quote}

\section{Implementation}
\label{sec:implementation}

The implementation of the algorithm is presented in two parts.  The
predictor step uses a suitable modification of the linear complementarity 
problem solver contained in the PATH code.  The corrector uses a heuristic
based on the Moore-Penrose idea.  We assume hereafter that $C = [\ell,u]$.

\subsection{Prediction}
\label{sec:pred-impl}
PATH \cite{dirkse.ferris:path} is an implementation of a nonsmooth Newton
method for solving mixed complementarity problems.  At each iteration,
a linear complementarity problem is solved using a homotopy method, i.e.
it follows the zero curve of a homotopy from $(x^k, 0)$ to $(x', 1)$.
The technique employed to generate this path is the same as that developed in
\secref{generating}.  
%The termination criteria is slightly different in 
%that $h_k$ is taken to be infinity, i.e. there is no restriction placed 
%upon the distance beween $x^k$ and $x'$.  
The code has 
special rules to deal with the cases where the Jacobian of the affine 
map on $\sigma$ does not have full row rank and where there is
degeneracy in the path \cite{ferris.munson:interfaces}.

The predictor step of our PL homotopy algorithm uses the same code to 
generate $Q^a_k$ for the homotopy $\LHC^a(\xk,\tk)(x,t)$ starting in a
cell $\sigma$.  
% The only difference
% is that we start from a point 
% $(x^k, t^k)$ where $\LHC^a(\xk,\tk)(\xk,\tk) \approx 0$.
We explain the implementation of the algorithm of \secref{predictor}
by generating a representation of this linearized
normal map, an indexing scheme that determines $\sigma$, 
and the Jacobian, $J$, of
the affine map restricted to $\sigma$.  We then discuss how a set of columns
spanning $J$ and the orientation is calculated.  
We continue with a presentation
of how the cells of the normal manifold are traversed and the termination
rules used by the code.

\subsubsection{Representation of the linearized normal map}
We will describe the linearized normal map 
$\LHC^a(x^k,t^k)(\cdot,\cdot\cdot)$ by giving a formula for the 
affine mapping that defines it on each cell.
We first decompose $x\in\real^n$ into a triple $(z,w,v) \in \real^{3n}$,
then give a simple indexing scheme for representing each cell of $\calN_C$,
and finally combine these two to describe $\LHC^a(x^k,t^k)(\cdot, \cdot\cdot)$
on each cell. 

Represent each $x$ as $z - w +v$ where
\begin{eqnarray}
\label{eqn:zwv}
z = \pi_C(x), \; w = [\pi_C(x) - x]_+, \; 
v = [x - \pi_C(x)]_+
\end{eqnarray}
and, for any vector $c$, $c_+$ is the vector whose $i$th component 
is $\max\{c_i,0\}$.  So $z\in C$ and $w$, $v$ satisfy 
\begin{eqnarray}
\label{eqn:wv ge 0}
w,v \ge 0
\end{eqnarray}
and the complementarity conditions
\begin{eqnarray}
\label{eqn:compl wv} 
\ip{z - l}{w} = 0, \;
\ip{u - z}{v} = 0.
\end{eqnarray}
Conversely, for any $z \in C$ and nonnegative 
$w$ and $v$ satisfying (\ref{eqn:compl wv}),
the vector given by $x = z-w+v$ is such that (\ref{eqn:zwv}) holds.

Define a {\em representation} $\calA = \{ \calA_i \}$ by
$\calA_i = i$, $n+i$ or $n+2i$, corresponding to $z_i$, $w_i$ or $v_i$
respectively, where 
$\calA_i \neq n+i$ if $l_i = -\infty$
and 
$\calA_i \neq n+2i$ if $u_i = \infty$.
We abuse notation by writing, say,
$w_i \in \calA$ to indicate that $n+i \in \calA$.
It is clear that there is a one-to-one correspondence between the
cells 
\[\begin{array}{rcccl}
\sigma_{\calA} &=& \{x \,: & x_i \le l_i & \mbox{if $w_i \in \calA$}
\\     & &         & l_i \le x_i \le u_i & \mbox{if $z_i \in \calA$}
\\     & &         & u_i \le x_i & \mbox{if $v_i \in \calA$} 
\}
\end{array}\] 
and the representations $\calA$, and that
$\calN_C$ is precisely the family of cells
indexed by representations $\calA$.

Let $\sigma_{\calA}$ be a cell and $\calA$ be its representation.
Then for $x \in \sigma_{\calA}$
and $(z,w,v)$ given by $(\ref{eqn:zwv})$, the relationships
(\ref{eqn:compl wv}) and (\ref{eqn:wv ge 0}) imply
\[ \begin{array}{lllllllllll}
z_i = l_i &\mbox{and}& w_i = l_i -x_i &\mbox{and} & v_i=0 & \mbox{if $w_i \in \calA$},
\\
z_i = x_i &\mbox{and}& w_i= 0 &\mbox{and}& v_i=0  &\mbox{if $z_i \in \calA$,}
\\
z_i = u_i &\mbox{and}& w_i=0 &\mbox{and} & v_i = x_i - u_i & \mbox{if $v_i \in \calA$}
.
\end{array}\]
Thus as $x$ varies within $\sigma_{\calA}$,
only the components of $(z,w,v)$ in the representation $\calA$ 
can change and, furthermore, 
the transformation $x \mapsto (z,w,v)$ is affine.

Given $x^k$ let $z^k = \pic{x^k}$, $w^k = [\pi_C(x^k) - x^k]_+$ 
and $v^k = [x^k - \pi_C(x^k)]_+.$
Then rewrite 
$\LHC^a(x^k,t^k)(x,t)=0$ as 
the equivalent system consisting of (\ref{eqn:wv ge 0}),
(\ref{eqn:compl wv}) together with
\begin{eqnarray}
\label{eqn:linear zwv}
N(z,w,v,t) + q & = & 0 ,
\end{eqnarray}
where $I \in \real^{n\times n}$ is the identity,
\begin{eqnarray*}
\begin{array}{rclcl}
N & = & [A \; -\!I \; I \; r] &\in& \real^{(3n+1)\times n}\\[3pt]
A & = & t^k\nabla f(z^k) + (1 - t^k)I  &\in& \real^{n\times n}\\[3pt]
r & = & f(z^k) - (z^k - a)  &\in& \real^n \\[3pt]
q & = & t^kf(z^k) + (1-t^k)(z^k - a) - Mz^k - t^kr  \\
& = & t^kz^k - a - t^k\nabla f(z^k)z^k &\in& \real^n .
\end{array}
\end{eqnarray*}
So a solution $(z,w,v,t)$ of (\ref{eqn:wv ge 0})--(\ref{eqn:linear zwv})
yields $\LHC^a(x^k,t^k)(z-w+v,t)=0$.

Since only the components in $\calA$ change when $x \in
\sigma_{\calA}$, it is clear 
that for all $x\in \sigma_{\calA}$ and $(z,w,v)$ representing $x$
\[ \LHC^a(x^k,t^k)(x,t) \;=\; N_{\calA}\,(z,w,v)_\calA \,+\, rt \,+\, q^\sigma.
\]
Here $q^\sigma = q + N_{\bar{\calA}}\, (z,w,v)_{\bar{\calA}}$
where 
$\bar{\calA} = \{1, \ldots, 3n\} \setminus \calA$
so that the subvector $(z,w,v)_{\bar{\calA}}$ is constant in $\sigma$.

The PATH code works with the triplet $(z,w,v)$ and $N_{\calA}$.
However, in $x$-space the Jacobian, $J$, of $\LHC^a(x^k,t^k)(x,t)$
restricted to $\sigma_{\calA}$ is given by
$\left[\begin{array}{cc} J_x & J_t\end{array}\right] =
\left[\begin{array}{cc} N_\calA\Delta & r\end{array}\right]$
where
$\Delta \in \Re^{n \times n}$ is the diagonal matrix
\[
\Delta_{i,i} =
\begin{cases}
-1 & \mbox{if $w_i \in \calA$} \\
1  & \mbox{otherwise.}
\end{cases}
\]
This fact will 
be used in the calculation of the orientation.

\subsubsection{Bases}
At the start of the predictor step at iteration $k$ we are given 
$x^k \in \real^k$, and a representation $\calA$ whose corresponding cell 
$\sigma$ contains $x^k$.  If $x^k$ is in the intersection of two or more
cells we use $\calA$ to describe the cell we are interested in.
We have seen that the linear part of $\LHC^a(x^k,t^k)$ on $\sigma$ can be 
described using the matrix 
$\left[\begin{array}{cc} N_\calA & r\end{array}\right] \in
\real^{n\times(n+1)}$.
To determine our initial search direction, we need to construct a
basis for the range space and kernel of this map.
More formally, given a cell $\sigma$ and its representation $\calA$, define a
{\em basis} $\calB$ as a list of $n$ indices taken from 
$\calA \cup \{3n+1\}$, 
where $3n+1$ is the index of $t$, such that $N_\calB$ is invertible.

The linear complementarity code is initially given a candidate basis,
$\tilde{\calB}$.  It is possible that 
$N_{\tilde{\calB}}$ is
not invertible.  In this case we attempt to uncover a
basis.  We apply a LU factorization to the $n \times (n+1)$ matrix
$N_{\calA \cup \{3n+1\}}$ to 
determine the linearly dependent column.  The remaining columns 
form a basis, $\calB$.

If more that one linearly dependent column is identified, the theory
breaks down and this procedure fails.  If $x^k$ is in
the intersection of two or more cells, we try to choose an adjacent cell 
with full row rank.  Specifically, the adjacent cell attempted is 
$\bar{\sigma}$ with representation:
\begin{eqnarray*}
\bar{\calA}_i =
\begin{cases}
w_i & \mbox{if $z^k_i = l_i$} \\
v_i & \mbox{if $z^k_i = u_i$} \\
z_i & \mbox{otherwise.}
\end{cases}
\end{eqnarray*}
We choose a candidate basis as $\bar{\calB} = \bar{\cal{A}}$. 
We then apply the above procedure with $\bar{\calB}$ replacing
$\tilde{\calB}$.  In the event this fails, an error is reported and
the code halts.

The reason we use a candidate basis as our initial guess as to the invertible 
basis as opposed to always factoring the larger $n \times (n+1)$ matrix is 
that the code written for PATH only uses the $n \times n$ basis
matrix for computations.  We did not want to rewrite the code to use
the $n \times (n+1)$ matrix.
Typically the candidate basis is invertible and we save a 
factorization.

\subsubsection{Orientation}
\label{sec:impl-orient}
We now have a initial basis $\calB^k$ and a linearly dependent column
$\beta^k = (\calA \cup \{3n + 1\}) \setminus \calB^k$.  We need to 
determine the initial direction of the path $Q^a_k$, that is whether the 
variable $(z,w,v,t)_{\beta^k}$ is going to initially increase or decrease. 
To do this, we need to find the orientation.  

For the orientation, we determine the sign of the determinant of the 
following matrix:
\begin{eqnarray*}
\begin{bmatrix}
J_x & J_t \\
d_{x}^T & d_{t}
\end{bmatrix} =
\begin{bmatrix}
N_\calA\Delta & r \\
d_x^T\Delta & d_t
\end{bmatrix} =
\begin{bmatrix}
B & b \\
d_B^T & 1
\end{bmatrix} P
\end{eqnarray*}
where $\left[\begin{array}{cc} d_x^T & d_t\end{array}\right]$ spans the kernel 
of $\left[\begin{array}{cc} J_x & J_t\end{array}\right]$, 
$B = \left[\begin{array}{cc} N_\calA\Delta & r \end{array}\right]_{\calB^k}$, 
$b = \left[\begin{array}{cc} N_\calA\Delta & r \end{array}\right]_{\beta^k}$, and
$P$ is a permutation matrix.  Let $LU$ be the decomposition of $B$ into
an invertible lower triangular matrix, $L$, and an invertible
upper triangular matrix, $U$.  Letting $d_B^T =  -U^{-1}L^{-1}b = -B^{-1}b$, 
we can see that $\left[\begin{array}{cc} d_B^T & 1 \end{array}\right]$
spans the kernel of $\left[\begin{array}{cc} B & b \end{array}\right]$.
Combining these facts, we have
\begin{eqnarray*}
\begin{bmatrix}
B & b \\
d_B^T & 1
\end{bmatrix}P
& = &
\begin{bmatrix}
LU & b \\
d_B^T & 1
\end{bmatrix}P \\
& = &
\begin{bmatrix}
L & 0 \\
d_B^TU^{-1} & 1 - d_B^TU^{-1}L^{-1}b
\end{bmatrix}
\begin{bmatrix}
U & L^{-1}b \\
0 & 1
\end{bmatrix}P \\
& = &
\begin{bmatrix}
L & 0 \\
d_B^TU^{-1} & 1 + \norm{d_B}^2
\end{bmatrix}
\begin{bmatrix}
U & L^{-1}b \\
0 & 1
\end{bmatrix}P.
\end{eqnarray*}
Therefore, the orientation is just the sign of the determinant of $L$ times 
the sign of the determinant of $U$ multiplied by the sign of determinant of $
P$.  We will initially
increase the entering variable if the orientation is $\eta$ and decrease it 
otherwise.  We then traverse the cells using a complementary pivoting 
strategy.

\subsubsection{Complementary Pivoting}
Given the basis, $\calB^k$, and the entering variable, $\beta^k$, we 
calculate a direction, 
$d = -N_{\calB^k}^{-1}
      N_{\beta^k}$.  We 
will either increase or decrease variable $(z,w,v,t)_{\beta^k}$ along this 
direction depending upon the orientation calculation above.  
A standard ratio test (on $(z,w,v,t)$) determines when we run into
the boundary of a cell and hence gives a
leaving variable and steplength.  We then perform a pivot by updating 
$\calB^k$ to include the entering variable and exclude the leaving variable.  
The new entering variable is chosen as follows:
\begin{itemize}
\item If $z_i$ leaves at its lower bound, $w_i$ is the new entering variable.
\item If $z_i$ leaves at its upper bound, $v_i$ is the new entering variable.
\item If $w_i$ leaves at its lower bound, $z_i$ is the new entering variable 
at its lower bound.
\item If $v_i$ leaves at its lower bound, $z_i$ is the new entering variable 
at its upper bound.
\item If $t$ leaves, we are done.
\end{itemize}
This defines the new cell along with a corresponding basis, entering
variable, and direction.  The code uses a rank-1 update of the LU 
factorization of the old basis to find the needed decomposition of the
new basis.  

\subsubsection{Termination}
The four conditions for termination presented in \secref{generating}
are used in the code.  If $t$ leaves the basis we terminate according to the 
complementary pivoting rules.  This corresponds to $t = 0$ or $t = 1$.  
In the code, we relax the distance 
based termination rule to be
\[ \norm{(\pi_C(x),t) - (\pi_C(x^k),t^k)} \;=\; h_k
\]
because the approximation error, i.e.~the difference between
$\LHC^a(x^k,t^k)(x,t)$ and $\HC^a(x,t)$
depends more on $\pic{x} - \pic{x^k}$
than on $[x-\pic{x}] - [x^k - \pic{x^k}]$.

\subsubsection{Special Rules}
The corrector can lead us to a place where $t^k > 1$.  In this case we reverse 
the direction we are traveling.  We put a lower bound on the $t$ variable of 
1 and an upper bound on $t$ of $\infty$.  We then use the above direction 
calculation, but take the opposite direction in which to initially travel.  
Everything then continues as normal, with the termination being that $t$ 
leaves at its lower bound or we take a step of the maximum allowed distance.

\subsection{Correction}
The corrector code is given an initial $(x',t')$ from the 
predictor.  The implementation of the corrector uses spacer steps which move 
the current iterate closer to $P^a$. They are defined as follows:
\begin{eqnarray*}
\hat{y}^0_i = \begin{cases}
(\pi_C(x') - [H^a(x',t')]_+)_i & \mbox{if variable $\pi_C(x') = l_i$}\\
(\pi_C(x') + [-H^a(x', t')]_+)_i & \mbox{if variable $\pi_C(x') = u_i$}\\
x_i' & \mbox{otherwise.}
\end{cases}
\end{eqnarray*}
Let $s^0 = t'$. 
We then calculate $\norm{H^a_C(y^0,s^0)}$ and check to see if it is zero
as detailed in \secref{heuristic}.  If it is not, we use the Moore-Penrose 
idea to move back onto the zero curve of the homotopy.

The Moore-Penrose iterate solves a minimization problem:
\begin{eqnarray*}
\begin{array}{ll} \mbox{min } & \norm{(x,t) - (y^j, s^j)}^2 \\
\mbox{subject to } &
\begin{bmatrix} J_x & J_t \end{bmatrix}
\begin{bmatrix} x \\ t\end{bmatrix} = -\HC^a(\hat{x},\hat{t})
\end{array}
\end{eqnarray*}
where $\left[\begin{array}{cc} J_x & J_t\end{array}\right]$ 
is the Jacobian of the affine map 
$\LHC^a(y^j, s^j)(\cdot, \cdot)$ on $\sigma$, a cell containing $y^j$,
(see \secref{pred-impl} for details on how this is calculated).

To solve this problem, we solve the following system of equations:
\begin{eqnarray*}
\begin{bmatrix}
I & 0 & -J_x^T \\
0 & 1 & -J_t^T \\
J_x & J_t & 0
\end{bmatrix}
\begin{bmatrix}
d_x \\
d_t \\
d_u
\end{bmatrix}
= \begin{bmatrix}
0 \\
-\HC^a(y^j,s^j)
\end{bmatrix}.
\end{eqnarray*}
The direction $(d_x, d_t)$ is nonzero since $\HC^a(y^j,s^j) \not = 0$.
We move in the direction $(d_x, d_t)$ until either the full step is
taken or we encounter the boundary of $\sigma$ to obtain $y^{j+1}$ and
$s^{j+1}$ and an adjacent cell $\sigma'$.  We the check for termination,
and if necessary perform another iteration.

The code does not allow the residual to increase from one corrector step to
the next and returns an error code if at the end the homotopy parameter
is negative. If the corrector fails, we take a half predictor step and
attempt the correct again.  A half predictor step is quickly found using
the reconstruction technique documented in \cite{ferris.munson:interfaces}.

\subsection{Updates}
The parameters in the code are updated after each major
iteration (predictor, corrector sequence).  We modify the maximum distance 
allowed to travel to become larger if we are doing well and make it smaller
if we are doing poorly.

\section{Results}
\label{sec:results2}

The algorithm was implemented using the current version of the linear
complementarity problem solver found in PATH.  The framework reported 
in \cite{ferris.munson:interfaces} was used so that we could easily 
access problems generated by the GAMS modeling language.  

We ran the code on the MCPLIB \cite{dirkse.ferris:mcplib} test suite of 
problems.  The MPSGE \cite{rutherford:applied} 
models in the collection were omitted due to the
fact that function and Jacobian evaluations for these models have the
undesirable side-effect of changing the evaluation point, thereby
making corrector steps impossible to perform.

The model name, dimension of the problem, number of
function evaluations (Func), and solution time (Time),
for the first starting point of the remaining models in this test set
is reported in Table \ref{mcplib}.  The test was carried out on a 
Sun Ultrasparc machine that has 768 megabytes of available memory.  

We note that most of the academic test problems which 
cause problems for other algorithms (e.g. billups, dirkse1, ralph,
simple-ex) are easily solved using the 
predictor-corrector code.
Over the entire set
of 77 models with 436 total starting points, we achieved
a 67.7\% success rate.  Of the 141 failures, 83 of these occurred for
starting points associated with just two problems, namely
\texttt{games} and \texttt{tinsmall}.  

\begin{table}
\footnotesize%\scriptsize
\begin{center}
\begin{tabular}{lrrr|lrrr}
Problem & Size & Func & Time & Problem & Size & Func & Time\\
\hline
badfree & 5 & 44 & 0.03 &       mr5mcf & 350 & - & fail\\
bert\_oc & 5000 & - & fail &    munson3 & 2 & - & fail\\
bertsekas & 15 & 441 & 0.66 &   munson4 & 1 & - & fail\\
billups & 1 & 32 & 0.02 &       nash & 10 & 203 & 0.71\\
bratu & 5625 & 181 & 13208.96 &         ne-hard & 3 & 111 & 0.08\\
cammcf & 242 & - & fail &       obstacle & 2500 & - & fail\\
choi & 13 & 21 & 2.66 &         opt\_cont & 288 & 173 & 12.10\\
colvdual & 20 & 269 & 0.45 &    opt\_cont127 & 4096 & 641 & 1726.96\\
colvnlp & 15 & 275 & 0.33 &     opt\_cont255 & 8192 & 1888 & 12470.15\\
cycle & 1 & 23 & 0.01 &         opt\_cont31 & 1024 & 323 & 114.73\\
degen & 2 & 27 & 0.02 &         opt\_cont511 & 16384 & - & fail\\
dirkse1 & 2 & 135 & 0.08 &      pgvon105 & 105 & 619 & 12.37\\
duopoly & 63 & - & fail &       pgvon106 & 106 & - & fail\\
eckstein & 1 & 21 & 0.01 &      pies & 42 & 1595 & 4.91\\
ehl\_k40 & 41 & 195 & 6.68 &    powell & 16 & 61 & 0.20\\
ehl\_k60 & 61 & 233 & 18.26 &   powell\_mcp & 8 & 75 & 0.11\\
ehl\_k80 & 81 & 265 & 38.50 &   qp & 4 & 39 & 0.03\\
ehl\_kost & 101 & 291 & 74.67 &         ralph & 2 & 33 & 0.02\\
electric & 158 & - & fail &     runge & 1 & 13 & 0.01\\
explcp & 16 & 98 & 0.19 &       scarfanum & 13 & 88 & 0.23\\
fixedpt & 406 & - & fail &      scarfasum & 14 & 46 & 0.13\\
forcebsm & 184 & - & fail &     scarfbnum & 39 & 653 & 2.46\\
forcedsa & 186 & - & fail &     scarfbsum & 40 & 621 & 3.28\\
freebert & 15 & 439 & 0.65 &    shubik & 30 & 1534 & 2.78\\
gafni & 5 & 48 & 0.07 &         simple-ex & 17 & 173 & 0.67\\
games & 16 & - & fail &         simple-red & 13 & 121 & 0.65\\
hanskoop & 14 & 58 & 0.10 &     sppe & 27 & 447 & 0.80\\
hansmcf & 43 & 139 & 0.74 &     tinloi & 146 & 54 & 1.16\\
hydroc06 & 29 & - & fail &      tinsmall & 42 & 647 & 3.58\\
hydroc20 & 99 & - & fail &      tobin & 42 & 305 & 0.79\\
jel & 6 & 317 & 0.37 &  trafelas & 2376 & - & fail\\
josephy & 4 & 43 & 0.03 &       trig & 10 & - & fail\\
kojshin & 4 & 61 & 0.04 &       vonthmcf & 125 & - & fail\\
kyh-scale & 28 & - & fail &     xu1 & 8 & 73 & 0.07\\
kyh & 28 & - & fail &   xu2 & 8 & 73 & 0.07\\
lincont & 419 & - & fail &      xu3 & 11 & 63 & 0.08\\
mathinum & 3 & 47 & 0.04 &      xu4 & 11 & 39 & 0.04\\
mathisum & 4 & 35 & 0.03 &      xu5 & 20 & 70 & 0.12\\
methan08 & 31 & - & fail
\end{tabular}
\caption{Results on MCPLIB Problems}\label{mcplib}
\end{center}
\normalsize
\end{table}

The current implementation demonstrates that homotopy methods can be
generated for large scale complementarity problems.  The choice of
homotopy function (in our case $x-a$) is critical for the success of
the these methods.  We were disappointed with the codes robustness; in
particular the implementation is very sensitive to how the original
problem is formulated (e.g. where $f_1(x) = 0$ or $-f_1(x) = 0$.  

Three approaches may be useful in improving these results.  Firstly, a
different choice of homotopy could be used instead of $x-a$ to
attempt to match it mored closely to the underlying problem.
Secondly,
the preprocessor for MCP described in
\cite{ferris.munson:preprocessing} could be used to exploit more fully
underlying problem structure.  Finally, this preprocessor could be
adapted to automatically identify poorly posed problems and
reformulate appropriately.  This is the subject of future research.

%\cite{ralph:global}
%\cite{robinson:mathematical}
%\cite{robinson:implicit-function}
%\cite{robinson:homeomorphism}
%\cite{robinson:nonsingularity}
%\cite{subramanian:gauss-newton}

%\bibliographystyle{plain}
%\bibliography{mathprog}

\medskip

{\noindent \large\bf Appendix}
\newcommand{\hatp}{\hat{p}}
\newcommand{\hatP}{\hat{P}}

We use a property of subanalytic sets \cite{lojasiewicz:ensembles,hironaka:introduction}
in the proof the next result:
If $S$ is a subanalytic set in $\real^n$
then for any point $\hatx$ in its closure there is a
continuous function $g:[0,1]\to \real^n$ with $g(1) = \hatx$
and $g(s) \in S$ for $s \in [0,1)$.

\begin{lemma}
If $p:[0,1) \to \real^n$ is a subanalytic function
then either there exists the limit $\lim_{t\to 1^-} p(t)$,
or $\norm{p(t)} \to \infty$ as $t \to 1^-$.
\end{lemma}

\begin{pf}
Let $p:[0,1) \to \real^n$ be subanalytic,
and $\hatp$ be a limit point of $p(t)$
as $t\to 1^-$. It is sufficient to show that 
the limit $\lim_{t\to 1^-} p(t)$ exists and equals~$\hatp$.

The graph of $p$, $P = \{(t,p(t))\,:\, t \in [0,1)\}$, 
is subanalytic by definition and $(1,\hatp)$ lies in its closure. 
Let %$\hatP = P \cup \{(1,\hatp)\}$, and 
$g:[0,1] \to \real^n$ be a continuous function with %$g(0) = (0,p(0))$ 
$g(1) = (1,\hatp)$ and $g(s) \in P$ for $s \in [0,1)$.
Denote the first component function
of $g(s)$ by $\gamma(s) \in [0,1]$.
Hence for $s \in [0,1)$, $\gamma(s) \in [0,1)$ and
$g(s) = (\gamma(s), p\circ\gamma(s))$.
%Define $\hats = \inf\{s \in [0,1] \,:\, \gamma(s) = 1 \}$.
%Then $\hats \in (0,1]$  
%and
%$\gamma(\hats) = 1$
%by continuity of $g$.
%Also  for each $s \in [0,\hats)$ we have
%$\gamma(s) \in [0,1)$ so that $g(s) \in P$, hence
%$g(s) = (\gamma(s), p\circ\gamma(s))$.
%Similarly, $\gamma(\hats) = 1$ implies $g(\hats) = (1,\hatp).$ 
% Define $\hats = \inf\{s \in [0,1] \,:\, g(s) = (1,\hatp) \}$.
% Then $\hats \in (0,1]$ 
% and
% $g(\hats) = (1,\hatp)$
% by continuity of $g$; 
% %$\gamma(\hats) = 1$;
% and for each $s \in [0,\hats)$ we have
% (i)~$\gamma(s) \in [0,1)$, and (ii)~$g(s) \in
% \hatP\setminus\{(1,\hatp)\} = P$ so that 
% $g(s) = (\gamma(s), p\circ\gamma(s))$.

Let $\epsilon > 0$. We will complete the proof by providing $\delta
\in (0,1]$ such that 
  %
\begin{equation}\label{eqn:cty of p} 
   \norm{p(t) - \hatp} \; < \; \epsilon 
   \quad\mbox{for all } t \in [1-\delta, 1).
\end{equation} 
  %
Continuity of $g$ at $s=1$ yields $\delta_g \in (0,1]$ such 
\[
\norm{p\circ\gamma(s) - \hatp} \; < \; \epsilon 
\quad\mbox{for all } s \in [1-\delta_s,1). 
\]
%
Now $\gamma(1-\delta_g) < 1=\gamma(1)$ so continuity
of $\gamma$ ensures (by the intermediate value theorem) that
$\gamma([1-\delta_g,1)) \supset [\gamma(1-\delta_g), 1)$.
Let $\delta = 1 - \gamma(1-\delta_g) \in (0,1]$ and deduce
(\ref{eqn:cty of p}) from the previous bound.  
\end{pf}

An easy corollary extends a result of 
\cite[Theorem~5.10]{facchinei.pang:total},
in the context of a homotopy method applied to
variational inequalities, %, {\em Tikhonov regularization},
in which parts~2 and~3 of were shown to be equivalent.
\begin{corollary} For a subanalytic function $p:[0,1) \to \real^n$
the following statements are equivalent:
\begin{enumerate}
\item
$\liminf_{t\to1^-} \norm{p(t)} < \infty$.
\item
$\limsup_{t\to1^-} \norm{p(t)} < \infty$.
\item
$\lim_{t\to1^-} p(t)$ exists.
\end{enumerate}
\end{corollary}
