\chapter{SEMI}\label{chap:semi}

PATH is currently the most widely used code for solving complementarity 
problems.  While it may be argued that piecewise linear maps are more 
effective at approximating piecewise smooth maps, generating the ``Newton'' 
step typically involves the arduous task of solving a linear complementarity 
problem.  A seemingly more attractive approach is to use an algorithm based 
on solving a system of linear equations to generate each ``Newton'' step.  
Recent theoretical work has outlined a host of methods with this property
\cite{
chen.harker:smooth,
chen.mangasarian:smoothing,
deluca.facchinei.ea:semismooth,
harker.xiao:newtons,
pang:b-differentiable,
pang.qi:nonsmooth}.
Amongst these, the semismooth algorithm \cite{deluca.facchinei.ea:semismooth} 
appears to have some of the strongest associated theory.

SEMI is the result of an effort to develop a code based upon the semismooth 
algorithm.  This chapter begins by briefly discussing the theoretical 
foundations of the semismooth algorithm in \secref{foundation}.  Many of the 
results contained in this section are given without proof; instead, we provide
references to the relevant literature. We then present the implementation 
details of the code.  The main focus is on the numerical aspects of the code 
used to overcome problems with singularity and ill-conditioning, and strategies
to recover from finding non-optimal stationary points of the merit function. 
Issues related to the use of iterative solvers for large scale problems 
are then outlined.  We demonstrate the code's robustness on the problems
in the GAMSLIB and MCPLIB \cite{dirkse.ferris:mcplib} collections by 
performing two tests, one using direct solution methods and another using 
only iterative techniques.  Both indicate that 
the SEMI code is reliable and scalable.

\section{Mathematical Foundation}\label{sec:foundation}
The semismooth solver is based on a reformulation of the mixed complementarity
problem as a nonlinear system of equations.  We first recall that a mapping
$ \phi : \Re^2 \to \Re $ is called an {\em NCP-function} if it satisfies
\[
   \phi (a,b) = 0 \Longleftrightarrow 0 \leq a \;\perp\; b \ge 0.
\]
Two examples of NCP-functions are the Fischer-Burmeister 
\cite{fischer:special} function
\begin{equation}\label{FB}
   \phi_{FB} (a,b) := \sqrt{a^2+b^2}-a-b
\end{equation}
and the penalized Fischer-Burmeister \cite{chen.chen.ea:penalized} function
\begin{equation}\label{CCK}
   \phi_{CCK} (a,b) := \lambda \left( \sqrt{a^2+b^2}-a-b \right) -
   (1- \lambda) \max \{ 0, a \} \max \{ 0, b \},
\end{equation}
where $ \lambda \in (0,1) $ is a given parameter.  To reformulate the mixed 
complementarity problem using these functions 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.  Hence, the subscripts in the 
above index sets indicate which bounds are finite, with the only exception 
of $ I_f $ which contains the free variables.

Let $ \phi_1$ and $ \phi_2 $ belong to $\{\phi_{FB}, \phi_{CCK}\}$.
Then we can extend the idea in \cite{billups:algorithms} and define an 
operator $ \Phi : \rn \to \rn $ componentwise as follows:
$$
   \Phi_i (x) := \left\{ \begin{array}{ll}
   \phi_1 (x_i- l_i, F_i (x)) & \mbox{if } i \in I_l, \\
   - \phi_1 (u_i-x_i, - F_i (x)) & \mbox{if } i \in I_u, \\
   \phi_2 (x_i - l_i, \phi_1 (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}
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}
We remark that \chapref{path} used the same reformulation of the 
complementarity problem for globalization purposes, but with 
$\phi_1 = \phi_2 = \phi_{FB}$.  The above definition of $\Phi$ is more
general in that we can choose any combination of $\phi_{FB}$ and 
$\phi_{CCK}$ for $\phi_1$ and $\phi_2$.

Note that $\Phi$ is not differentiable in general. However, if $F$ is
continuously differentiable then $\Phi$ is semismooth. A standard technique to 
solve the mixed complementarity problem is to apply a nonsmooth Newton method 
(see \cite{qi.sun:nonsmooth,qi:convergence}) to the system $ \Phi (x) = 0 $ 
and globalize it using the corresponding merit function
\[
   \Psi (x) := \frac{1}{2} \Phi (x)\trt \Phi (x) = \frac{1}{2}
   \| \Phi (x) \|^2.
\]
Assuming that $ \Psi $ is continuously differentiable and recalling that the 
{\em B-subdifferential} of $ \Phi $ at a point $ x \in \rn $ is defined by 
\cite{qi:convergence}
\[
   \partial_B \Phi (x) := \{ H \in \rntn \, | \, \exists \{ x^k \}
   \subseteq D_{\Phi}: x^k \to x \mbox{ and } \Phi'(x^k) \to H \},
\]
where $ D_{\Phi} $ denotes the set of differentiable points of $ \Phi $,
we can follow the pattern from \cite{deluca.facchinei.ea:semismooth} and write 
down the basic semismooth solver for MCP.

{\bf Algorithm BSM (Basic Semismooth Method)}
{\it
\begin{itemize}
   \item[(S.0)] (Initialization) \\
      Choose $ x^0 \in \rn, \rho > 0, \beta \in (0,1), \sigma \in (0,1/2),
      p > 2, $ and set $ k:=0. $
   \item[(S.1)] (Stopping Criterion) \\
      If $ x^k $ satisfies a suitable termination criterion: STOP.
   \item[(S.2)] (Search Direction Calculation) \\
      Select an element $ H_k \in \partial_B \Phi(x^k). $ Find a solution
      $ d^k \in \rn $ of the linear system
      \begin{equation}\label{NEN}
         H_k d = - \Phi(x^k).
      \end{equation}
      If this system is not solvable or if the descent condition
      \begin{equation}\label{DCN}
         \nabla \Psi (x^k)\trt d^k \le - \rho \| d^k \|^p
      \end{equation}
      is not satisfied, set $ d^k:= - \nabla \Psi (x^k). $
   \item[(S.3)] (Line Search) \\
      Compute $ t_k:= \max \{ \beta^i \, | \, i = 0, 1, 2, \ldots \} $
      such that
      $$
         \Psi (x^k+t_k d^k) \le \Psi (x^k)+ t_k \nabla \Psi(x^k)\trt d^k.
      $$
   \item[(S.4)] (Update) \\
      Set $ x^{k+1}:=x^k+t_k d^k, k \leftarrow k+1, $ and go to (S.1). 
\end{itemize}
}

\noindent
Note that Algorithm BSM actually represents a whole class of methods
since it depends on the definition of $ \Phi $ which, in turn,
is completely determined by the choice of
$ \phi_1 $ and $ \phi_2 $. Note that, usually, $ \phi_1 $ plays the
central role in the definition of $ \Phi $; for example,  if there is no
variable with both finite lower and upper bounds, then $ \phi_2 $
is not used. In particular, this is the
case for the standard nonlinear complementarity problem.

For the purpose of this chapter, we are particularly interested in the
following two choices of $ \Phi $.  We define
$$
   \Phi_{FB} := \Phi \quad \mbox{if} \quad \phi_1 = \phi_2 = \phi_{FB},
$$
and
$$
   \Phi_{CCK} := \Phi \quad \mbox{if} \quad \phi_1 = \phi_{CCK}, \, \phi_2 =
   \phi_{FB}.
$$
The reader may wonder why we do not take $ \phi_2 = \phi_{CCK} $ as well; the 
simple reason is that we were unable to prove some of the subsequent results
for this case. In fact, basically all of these results are based on a suitable
overestimation for the generalized Jacobian $ \partial \Phi (x) $. Typically, 
such an overestimation can be obtained by exploiting Theorem 2.3.9 in 
\cite{clarke:optimization} that contains a convex hull operation. 
It is often possible to remove this convex hull and to get a simpler 
overestimate for $ \partial \Phi (x) $. However, when using 
$ \phi_1 = \phi_2 = \phi_{CCK} $, we were not able to remove the convex hull, 
so we decided not to take $ \phi_2 = \phi_{CCK} $ in the definition of the 
operator $ \Phi_{CCK} $.

Both from a theoretical and a numerical point of view 
\cite{chen.chen.ea:penalized}, the operator $ \Phi_{CCK} $ has stronger 
properties than $ \Phi_{FB} $, at least for the standard nonlinear 
complementarity problem.  Hence, $ \Phi_{CCK} $ will be used by default in the 
implementation of Algorithm BSM. However, in some situations, it is also 
helpful to have alternative operators like $ \Phi_{FB} $. For example,  
the implementation uses $\Phi_{FB}$ in one of the restarts.

We now summarize some of the properties of
$ \Phi_{FB} $ and $ \Phi_{CCK} $ as well as of their corresponding
merit functions
$$
   \Psi_{FB}(x):= \frac{1}{2} \Phi_{FB}(x)\trt \Phi_{FB}(x) \quad
   \mbox{and} \quad
   \Psi_{CCK}(x):= \frac{1}{2} \Phi_{CCK}(x)\trt \Phi_{CCK}(x).
$$
The proofs of the results can be found in \cite{billups:algorithms} and
\secref{theory} for the case of $ \Phi = \Phi_{FB} $. Since the proofs for 
$ \Phi = \Phi_{CCK} $ are very similar (although quite technical and lengthy),
we skip the proofs of all these results here.

\begin{proposition}\label{Semismooth}
Let $ \Phi$ belong to $ \{ \Phi_{FB}, \Phi_{CCK} \} $ and $\Psi$ be the 
corresponding merit function. Then the following hold:
\begin{itemize}
\item[(a)] $ \Phi $ is semismooth.
\item[(b)] If $ F' $ is locally Lipschitzian, then $ \Phi $ is strongly semismooth.
\item[(c)] $ \Psi $ is continuously differentiable on $ \rn $.
\item[(d)] If $ x^* $ is a strongly regular solution of MCP, then 
$ x^* $ is a BD-regular solution of $ \Phi (x) = 0 $.
\end{itemize}
\end{proposition}

\noindent
We only note that (strong) semismoothness is one of the two ingredients which 
are needed to prove local superlinear (quadratic) convergence of a nonsmooth 
Newton method. The second ingredient is the BD-regularity assumption.
See \chapref{introduction} for precise definitions of these terms.

We can now state the following convergence properties 
of Algorithm BSM. The proof is analogous to those given in 
\cite{deluca.facchinei.ea:semismooth} for $ \Phi_{FB} $ and the standard 
nonlinear complementarity problem.

\begin{theorem}\label{Convergence}
Let $ \Phi \in \{ \Phi_{FB}, \Phi_{CCK} \} $ and $ \{ x^k \} $ be a sequence 
generated by Algorithm BSM.  Then any accumulation point of this sequence is 
a stationary point of $ \Psi $. Moreover, if one of these accumulation points,
say $ x^* $, is a BD-regular solution of $ \Phi (x) = 0 $, then the
following statements hold:
\begin{itemize}
   \item[(a)] The entire sequence $ \{ x^k \} $ converges to $ x^*. $
   \item[(b)] The search direction $ d^k $ is eventually given by
      the Newton equation (\ref{NEN}).
   \item[(c)] The full stepsize $ t_k = 1 $ is eventually accepted in
      Step (S.3).
   \item[(d)] The rate of convergence is Q-superlinear.
   \item[(e)] If $ F' $ is locally Lipschitzian, then the rate of convergence 
              is Q-quadratic.
\end{itemize}
\end{theorem}

\section{The Linear System}
The heart of the semismooth algorithm lies in the linear algebra.  Most of
the time is spent solving the Newton system, $H_k d = -\Phi(x^k)$, where
in general $H_k$ is neither symmetric nor positive definite.  Effective 
mechanisms for solving this system using either iterative techniques or a 
direct method are indispensable and have great impact upon the success of 
the algorithm.  The key advantage of the semismooth algorithm over 
PATH 4.x  is that the former only solves a single linear 
system per iteration while the latter uses a pivotal based code to solve a 
linear complementarity problem.  The pivotal based code relies upon the 
availability of a direct factorization and rank-1 updates which limits its 
applicability to medium sized or large, structured problems.  Semismooth has
no such restriction.

We begin the analysis by investigating iterative techniques for finding the 
Newton direction as these will enable the algorithm to solve very large 
problems.  We present three of the methods considered and discuss time and 
space requirements and scaling issues.  The methods were evaluated by 
applying them to two reasonably large models representative of the test suite.
From the results, we conclude that LSQR \cite{paige.saunders:lsqr} is the most 
reliable.  Practical termination criteria are also mentioned.

We then discuss the issues involved and options available when using
direct methods.  The main difficulty encountered is singularity in the Newton 
system.  Information on detecting singularity and using that knowledge to 
construct a useful direction even in this case is presented.  The effects of 
the techniques considered on the singular models in the test set are given 
and the final choice of strategies for solving the linear system is provided.

Wherever possible, we want to use the best available technique to determine 
the Newton direction.  While LSQR is very reliable, typically the best 
technique in terms of time is to use a direct method to factor $H_k$.
However, when the size of the factors grows too large, we want to resort to 
the iterative technique.  Therefore, the rules implemented are designed so 
that large, structured problems with sparse decompositions will use the 
factorization software, while problems that are either very dense or have
large factors will not.  Currently, a restriction on the size of the
decomposition of 12 million nonzeros is imposed.  We no longer consider
using a direct method if this condition is violated.
%If the size of the decomposition contains more than 
%$b$ million nonzeros, we will no longer consider the direct
%factorization;
%currently $b$ is equal to  12 million.

\subsection{Iterative Techniques}
Three iterative techniques for finding the Newton direction were investigated:
LSQR, GMRES, and QMR.  We recall that the systems of equations solved will 
generally be neither symmetric nor positive definite.  Therefore, we cannot 
directly use the popular techniques from optimization algorithms such as 
conjugate gradients.  The algorithms tested are a representative set of those 
meeting the requirements.

LSQR \cite{paige.saunders:lsqr} is based upon the bidiagonalization method 
developed in \cite{golub.kahan:calculating} which implicitly solves the least 
squares problem, $\min{\norm{H_k d + \Phi(x^k)}^2}$.  The method is 
essentially a reliable variant of conjugate gradients applied to the normal 
equation, $H_k^TH_k d = -H_k^T\Phi(x^k)$.  The code only requires a workspace 
of 3 $n$-vectors in addition to the storage of $H_k$, $d$, and 
$\Phi(x^k)$.  The cost per iteration consists of $4nnz + 12n$ floating point 
operations, where $nnz$ is the number of nonzero elements in $H_k$.

The GMRES \cite{saad:iterative,walker:implementation} method uses the Arnoldi 
procedure to construct an orthonormal basis for the Krylov subspace
$\mathcal{K}_m\left(H_k, -\frac{\Phi(x^k)}{\norm{\Phi(x^k)}}\right)$, where
$$
   \mathcal{K}_m (A,r):= \mbox{span} \left\{ r, Ar, \ldots, A^{m-1}r
   \right\}
$$
for some matrix $ A \in \rntn $ and a vector $ r \in \rn $.
We then 
use this basis to find a vector in the generated subspace minimizing the 
residual, $\norm{H_k d + \Phi(x^k)}^2$.  The implementation uses Householder 
reflections
for the orthogonalization process to preserve stability, maintains the current
optimal value of the residual at each iteration using plane rotations, and 
restarts after 
$m$ iterations.  We remark that because the matrices are not guaranteed to be 
positive definite, restarted GMRES can stagnate and make no progress.
The method has a workspace requirement of ($m+2$) n-  and 4 m-vectors and 
uses $2nnz + 4n(1+2i) - 4i^2$ operations per iteration where $i \in [1, m]$ 
is the iteration number.  Furthermore, we use around $4mn$ operations at the 
end of each $m$ iterations to generate the minimizing vector.  The main 
difficulty with GMRES lies in choosing the restart frequency.   If it is too 
small, we can fail to converge entirely, and if it is too large, the per 
iteration cost and storage requirements become significant.

The QMR \cite{saad:iterative} algorithm uses the Lanczos 
biorthogonalization algorithm to 
construct bases for the Krylov subspaces
$\mathcal{K}_m\left(H_k, -\frac{\Phi(x^k)}{\norm{\Phi(x^k)}}\right)$ and 
$\mathcal{K}_m\left(H_k^T, -\frac{\Phi(x^k)}{\norm{\Phi(x^k)}}\right)$
satisfying 
a biorthogonality condition.  The bases generated are then used to find
a vector with approximate minimum residual in 
$\mathcal{K}_m\left(H_k, -\frac{\Phi(x^k)}{\norm{\Phi(x^k)}}\right)$.  The 
QMRPACK \cite{freund.nachtigal:qmrpack} code tested uses the coupled two-term 
recurrence variant of the look-ahead 
Lanczos algorithm.  We allowed $m$ look-ahead steps to be tried, which
results in a workspace of $10m + 1$ n- and $8m + 18$ m-vectors.  Typically, 
$m$ is chosen to be small.  The code uses at least $4nnz + 14n$ floating point 
operations per iteration and requires a backsolve at the end to determine the 
vector with approximate minimal residual.

Scaling the linear system is crucial to the success of the iterative 
algorithms.  We define a matrix scaling using the following diagonal matrices
$R,C\in\Re^{n\times n}$ with diagonal entries:
\begin{eqnarray}
\begin{array}{lcl}
R_{i,i} & = &
\frac{1}{\max\left\{\sqrt{\Phi_i(x^k)^2 + \sum_j (H_k)_{i,j}^2},
10^{-10}\right\}},\\
C_{j,j} & = & \frac{1}{\max\left\{\sqrt{\sum_i (RH_k)_{i,j}^2}, 
10^{-10}\right\}}.
\end{array}\label{row-scale}
\end{eqnarray}
We solve the linear system $RH_kCC^{-1}d = -R\Phi(x^k)$ by defining
$\tilde{d}:=C^{-1}d$, solving the system $RH_kC\tilde{d} = -R\Phi(x^k)$,
and recovering the Newton direction as $d = C\tilde{d}$.  
This procedure scales the rows and then the columns so that each has a two 
norm of 1.  The constants in the $\max$ operator are used to avoid division 
by zero errors.  When a row or column has a two norm close to zero, the 
scaling has no effect.  Scaling significantly reduces the number of 
iterations required in all cases and thus the total time spent in the 
iterative solver.

We remark that the use of preconditioned versions of LSQR, GMRES, and
QMR could further reduce the number of iterations taken and time spent to 
solve the system of equations.  This is a topic of future research.

\subsubsection{Evaluation}
We selected two reasonably large models with known solutions from the test 
suite on which to evaluate the iterative techniques.  The {\tt uruguay} model
has 2,281 
rows and columns and 90,206 nonzeros and is interesting because a direct 
factorization incurs a large amount of fill-in.  {\tt opt\_cont255} is a 
structured problem with 8,192 rows and columns and 147,200 nonzeros.  For each 
complementarity model, 
we only solved the first linear system arising from 
(\ref{NEN}).  The termination criteria 
for these tests was based upon the relative residual, 
$\frac{\norm{H_0 d_i + \Phi(x^0)}}{\norm{\Phi(x^0)}}$.  The iterative methods
terminated
when the relative residual is less than $10^{-8}$.  In all cases, we chose 
an initial guess of $d_0:=0$.  We did not investigate other choices.

When using the GMRES method, we need to choose the restart frequency.
We varied the value of $m$ by choosing values between 10 and 200.  The results
for each value of $m$ tested are given in the accompanying tables.  
%The
%same code was used in each case, we only varied the restart frequency be
%using the {\tt gmres\_restart\_frequency} option which has an upper bound of
%200.  
For QMR, we need to determine the number of look-ahead steps allowed.
We placed an upper limit of 20 on these steps, although smaller values
would suffice.

All of the trials were run on the same machine using the same executable
so that we can make a valid comparison.  
%Only the {\tt inexact\_method} 
%and {\tt gmres\_restart\_frequency} 
%options were varied during the tests.  Valid
%choices for {\tt inexact\_method} are {\tt lsqr}, {\tt gmres}, and {\tt qmr}.
Tables \ref{eval1} and \ref{eval2} report the results on the two test
problems.  The total iterations and time, exit status, and relative residual
at the final point are given.  Based upon 
the available information, we conclude that LSQR is the best choice for this 
purpose.  This method is quite effective for the models we have generated.  
However, it may require a large number of iterations in order to converge.  
In these situations we are willing to sacrifice speed in exchange for 
reliability.  Note that results shown later demonstrate the robustness of 
LSQR on the entire test suite. The robustness and effectiveness of
LSQR is also in accordance with the results of 
\cite{deluca.facchinei.ea:theoretical}.

\begin{table}
\begin{center}
\caption{Iterative Method Results: {\tt uruguay}}\label{eval1}
\begin{tabular}{l|cccc}
Method & Iterations & Time & Status & Relative Error \\
\hline
LSQR      &    908 &  43 & Solved & 1.0e-8 \\
GMRES 10  &  1,766 &  47 & Solved & 1.0e-8 \\
GMRES 20  &  1,346 &  46 & Solved & 1.0e-8 \\
GMRES 50  &    799 &  45 & Solved & 9.9e-9 \\
GMRES 100 &    481 &  45 & Solved & 9.7e-9 \\
GMRES 200 &    419 &  54 & Solved & 9.9e-9 \\
QMR       &    466 &  41 & Solved & 6.5e-9
\end{tabular}
\end{center}
\end{table}

\begin{table}
\begin{center}
\caption{Iterative Method Results: {\tt opt\_cont255}}\label{eval2}
\begin{tabular}{l|cccc}
Method & Iterations & Time & Status & Relative Error \\
\hline
LSQR      &   2,848 &     90 & Solved          & 9.9e-9 \\
GMRES 10  & 100,000 &  2,114 & Iteration Limit & 8.7e-4 \\
GMRES 20  & 100,000 &  3,194 & Iteration Limit & 2.6e-4 \\
GMRES 50  & 100,000 &  6,217 & Iteration Limit & 6.0e-5 \\
GMRES 100 & 100,000 & 11,498 & Iteration Limit & 1.1e-5 \\
GMRES 200 & 100,000 & 22,136 & Iteration Limit & 5.4e-8 \\
QMR       & 100,000 &  5,288 & Iteration Limit & 1.6e-7
\end{tabular}
\end{center}
\end{table}

\subsubsection{Termination Rules}\label{termination}
While the termination rule given above is reasonable for evaluation,
we now return to the subject of practical termination rules.  The relative
residual calculated above is not applicable as a termination criterion unless
we know a priori that the  linear model has a solution.
This is an unreasonable assumption to make.  Therefore, the implementation of
LSQR uses the termination rules developed in \cite{paige.saunders:lsqr}.
They are to terminate if any of the following holds:
\begin{enumerate}
\item $\mbox{cond}(H_k) \ge \mbox{CONLIM}$
\item $\norm{r_i} \le \mbox{BTOL}\norm{\Phi(x^k)} + \mbox{ATOL}\norm{H_k}\norm{d_i}$
\item $\frac{\norm{H_k^Tr_i}}{\norm{H_k}\norm{r_i}} \le \mbox{ATOL}$
\end{enumerate}
where $r_i = -(H_kd_i + \Phi(x^k))$.
Justification of these rules and a demonstration of their effectiveness is
given in \cite{paige.saunders:lsqr}.  We note that LSQR builds up estimates of
$\norm{H_k}$ and $\mbox{cond}(H_k)$ by performing a small amount of additional
computation per iteration of the code.  The exact tolerances used are 
$ATOL = \epsilon^{\frac{2}{3}}$, $BTOL = \epsilon^{\frac{2}{3}}$,
and $CONLIM = \frac{1}{10\sqrt{\epsilon}}$ where $\epsilon$ is the machine
precision.  Furthermore, an iteration limit of $\min\{100000, 20n\}$ was 
used.  These tolerances force the iterative method to find a point close 
to the exact solution 
of the linear system if it exists.  We did not investigate using less stringent
termination criteria.

\subsection{Direct Methods}
We now look at the issues involved in using direct
methods to solve the Newton system.  The factorization software
needs to have routines to factor and solve, and should be able
to uncover singularity problems and make a good approximation of the
linearly dependent rows and columns in such cases. For reasonably sized problems
we use the LUSOL \cite{gill.murray.ea:maintaining} sparse factorization 
routines contained in the MINOS \cite{murtagh.saunders:minos} nonlinear 
programming solver to factor $H_k$ and solve for the Newton direction.
The authors of this package have investigated the effects of modifying 
tolerances in the factorization on general linear systems and have suggested 
defaults which we have adopted for all the results.

The major difficulty with the direction finding problem is dealing with
those instances where the Newton system does not have a solution.  These 
singularity problems frequently occur in real world applications.  However, 
the theoretical algorithm only provides a crude mechanism in this case, i.e. 
the use of a gradient step, while other approaches may be more effective.
Clearly, any practical implementation of the semismooth algorithm must 
include appropriate procedures to deal with singularity.

Following the success of scaling for the iterative techniques, we first
investigate the applicability of scaling in conjunction with direct methods 
to avoid ill-conditioned systems.  We then look at techniques to determine a 
useful direction when the model is singular, including using gradient steps, 
diagonal perturbations of $H_k$, and finding least squares solutions to the
linear system.  Empirical evidence is provided upon which we evaluate the
methods.

\subsubsection{Scaling}
LUSOL contains routines to detect when a matrix is singular or nearly 
singular.  We study in this subsection the effects of scaling the linear 
problems in an effort to improve the conditioning of the matrices that we 
request to factor.  We want to determine if we can significantly reduce the 
number of occurrences where the factorization package determines that the 
matrix is singular.  By using scaling we hope to improve the overall
reliability of the code on ill-conditioned problems.

Two different scaling schemes were tested on the problems in the test 
set along with the default of no scaling.  The first technique is the diagonal
scaling used in the PATH 4.x solver.  In this case, we define a row scaling by
looking at elements of the diagonal of $H_k$ which are large and scale the
entire row of the problem.  Formally, we define a diagonal matrix $R$ such
that if $|(H_k)_{i,i}| > 100$ then 
$R_{i,i} = \frac{10}{|(H_k)_{i,i}|}$ and $R_{i,i} = 1$ otherwise.  We 
then try to factorize the scaled matrix $RH_k$.  
%This type of scaling is 
%turned on by using the option {\tt scale\_method bigdiagonal}.

The other scaling method is the matrix scaling as used in the iterative
techniques.  We use the diagonal matrices defined in (\ref{row-scale})
and attempt to factorize $RH_kC$.  
%This type of scaling is turned on by using 
%the option {\tt scale\_method matrix}.

There are costs associated with scaling. Of the two methods, matrix scaling 
is more expensive per iteration because it requires looking at the data twice.
We tested all of these scalings on the models in the entire test set and 
report in Table \ref{scaling} the number of detected singular solves, failures
of the algorithm to find a solution, and total time in seconds over the
entire test set.  
%Only the {\tt scale\_method} option was changed during
%this test; all of the other options were held constant.
When a singular model was detected we use the least squares recovery method 
detailed in Section \ref{recovery}.  Since diagonal scaling does not improve
significantly over no scaling, we disregard this method. 
The reason for this poor behavior is probably due to the fact that the
diagonal elements do not necessarily reflect the actual scaling of the
problem.  Matrix scaling 
significantly reduces the number of singular systems detected.  However, it 
does result in additional failures of the algorithm.  
Since we were unable to definitively choose between no scaling and
matrix scaling,
in the next section we
look at recovery techniques and report results for both.

\begin{table}
\begin{center}
\caption{Scaling Effects on Direct Methods}\label{scaling}
\begin{tabular}{l|ccc}
 & Detected \\
Scaling & Singular Matrices & Failures & Time \\
\hline
None & 2138 & 28 & 13,922 \\
Diagonal & 2248 & 30 & 12,780 \\
Matrix & 1331 & 30 & 12,901 \\
\end{tabular}
\end{center}
\end{table}

Before continuing, we note that the scaling investigated here is not very 
exhaustive and more complex schemes might be tested.  Furthermore, the 
scaling is being performed on the linear model, when it might be more 
appropriate to look at the nonlinear model to determine the scaling.
Finally, we did not investigate modifying other parameters, such as those
encountered in the nonlinear model in Section \ref{nonlinear}, in 
conjunction with scaling, which might lead to improved reliability and 
performance.

\subsubsection{Singularity}\label{recovery}
Having looked at scaling we need to establish procedures to recover from the 
singularity problem and generate a reasonable direction.  We have 
investigated three techniques.  The first is the theoretical standby of
using only gradient steps when the Newton system is unsolvable.  A second
technique is to use a diagonal perturbation of $H_k$ to regularize the problem.
The final method is to use LSQR to calculate a least squares solution of the
system.  All of the results in this section are only given for those models 
where singularity was detected by the linear solver.  Only a few options
were changed for each run, with the rest being held constant.

\paragraph{Gradient Steps}
The naive recovery technique, and the simplest of those considered, is to
simply resort to a gradient step whenever a singular model is detected.
This approach is theoretically justified, but in practice frequently leads 
to a stationary point that does not solve the complementarity
problem.  
However, we use this approach as the baseline against 
which we evaluate the rest of the methods.  The results for the 
194 singular models are given in Table
\ref{gradient} where we report the number of times the algorithm failed and 
total time for both scaled and unscaled models.  In this case, we note that 
matrix scaling performs better than no scaling.

\begin{table}
\begin{center}
\caption{Gradient Results on Singular Models}\label{gradient}
\begin{tabular}{l|cc}
Scaling & Failures & Time \\
\hline
none   & 28 & 4,836 \\
matrix & 24 & 3,630
\end{tabular}
\end{center}
\end{table}

\paragraph{Perturbation}
Perturbation involves replacing the linear model with one which
does not have a singularity problem.  We investigated using a diagonal
perturbation where we replace $H_k$ with 
$H_k + \lambda I$ for some $\lambda > 0$.  The value for $\lambda$ was chosen 
in the interval $[\alpha, \beta]$, with $\lambda = \gamma\Psi$ whenever 
possible.  
We used values of $\alpha = 10^{-8}$, $\beta = 1$ and $\gamma = \frac{1}{10}$.  
%The values of these constants can be changed with the 
%{\tt newton\_perturbation\_minimun}, {\tt newton\_perturbation\_maximum}, and 
%{\tt newton\_perturbation\_merit\_factor} options respectively.
When the perturbation is insufficient to overcome the singularity,
we increase $\lambda$ to $\delta\lambda$ for some $\delta > 1$.
We currently use $\delta = 10$, and  allow the perturbation
to increase only one time per iteration.
% which can be modified with the
% {\tt newton\_perturbation\_increase\_factor} option.  The number of times
%we allow the perturbation to increase per iteration can be controlled by
%the {\tt newton\_perturbation\_maximum\_iterations} which has a default value
%of 1.

The other choice to make is when to add the perturbation.  There are two
options investigated:
\begin{itemize}
\item If the first model encountered is singular, calculate a $\lambda$ and
monotonically decrease it from one iteration to the next.  The new value 
is $\min\left\{\kappa_1\lambda, \kappa_2\Psi\right\}$ where 
$\kappa_1 = 0.4$ and
$\kappa_2 = 0.1$ by default.  This strategy is used in the 
PATH 4.x code.  
%An initial perturbation can be set with the 
% {\tt proximal\_perturbation} option, $\kappa_1$ can be set with
% {\tt proximal\_perturbation\_shrink\_factor}, and $\kappa_2$ with
% {\tt proximal\_perturbation\_merit\_factor}.
\item Every time a singular model is encountered, calculate a value for
$\lambda$.
\end{itemize}
%These methods are specified with the options 
% {\tt proximal\_perturbation\_first yes} and
% {\tt proximal\_perturbation\_demand yes} respectively.  Setting both of these
%options results in a hybrid of the two which was not investigated. 
When  scaling was used, we first perturbed the problem and then scaled it.  

To test these strategies, we ran the set of singular models using each of
the options.  We report the number of failures on the 194 singular models
and total time in Table \ref{perturb}.  When the perturbation fails to find a
nonsingular matrix, the least squares method (to be described in the
next section) was used to calculate a 
direction.
The use of perturbation leads to fewer total failures than only using
the gradient method.  The effects of scaling the problem are mixed, with 
scaling leading to a decrease in total time, but sometimes resulting in 
additional total failures.

\begin{table}
\begin{center}
\caption{Perturbation Effects on Singular Models}\label{perturb}
\begin{tabular}{ll|cc}
Scaling & Strategy & Failures & Time \\
\hline
none   & first  & 21 & 8,910 \\
       & demand & 21 & 3,285 \\
matrix & first  & 20 & 8,210 \\
       & demand & 23 & 2,703 
\end{tabular}
\end{center}
\end{table}

\paragraph{Least Squares Method}\label{least}
Finally, we investigate the use of the LSQR iterative scheme to find a
solution to the least squares problem $\min{\norm{H_kd + \Phi(x^k)}^2}$ and
use the resulting $d$ as the Newton direction.  The practical termination 
rules mentioned in \ref{termination} were used.  

We investigated using scaling and no scaling in the linear model that we try 
to factor.  We present the results on the 194 singular models in Table
\ref{iterative} where we report the total number of failures in the algorithm
and time.  The major downside to using the iterative technique to solve the
least squares problem is that it is fairly slow because we allow many
iterations and have low tolerances.  We did not study the effect of changing 
the termination criteria.  However, the results indicate that this method is 
better than the others tested in terms of reliability.

\begin{table}
\begin{center}
\caption{LSQR Results on Singular Models}\label{iterative}
\begin{tabular}{l|cc}
Scaling & Failures & Time \\
\hline
none   &  9 & 8,513 \\
matrix & 11 & 7,534 \\
\end{tabular}
\end{center}
\end{table}

\subsection{Summary}
The empirical results given above provide clear choices.  For both large 
scale work and calculating a direction when the Newton system is singular 
we will use the LSQR iterative technique.  We remark that while this is 
the most reliable choice, it is perhaps not the most efficient method.
While we always scale the linear system when an iterative solver is used,
the effects of scaling the matrix we try to factor are indeterminate and
we made the decision to use no scaling to achieve simplicity in the
code.  
The effect of choices made in the nonlinear model which are
discussed in the next section have a great impact upon the success of the 
algorithm.  However, we did not investigate modifying those strategies in 
conjunction with the strategies in the linear solver.

%We provide a comprehensive listing of the options affecting the linear
%system solves in Table \ref{compr}.  We list the option name, default value,
%and a brief description of the purpose.

%\begin{table}
%\begin{center}
%\begin{tabular}{lcp{2.5in}}
%Option & Default & Description \\
%\hline
%decomposition\_size\_limit & 12 & Maximum number of nonzeros allowed 
%in factorization\\
%gmres\_restart\_frequency & 20 & Orthonormal 
%vectors created before restarting in GMRES \\
%inexact\_maximum\_iterations & 100,000 & 
%Absolute maximum number of iterations allowed in iterative solver\\
%inexact\_method & lsqr & Iterative solver used \\
%inexact\_scale\_method & matrix & Scaling used in the iterative solver \\
%inexact\_size\_factor & 20 & Factor of problem size to use as
%iteration 
%limit in iterative solver\\
%newton\_perturbation\_increase\_factor & 1 & Fraction used to increase the
%perturbation when necessary \\
%newton\_perturbation\_maximum & 1 & Maximum perturbation allowed in 
%Newton system \\
%newton\_perturbation\_maximum\_iterations & 1 & Maximum iterations of 
%perturbation scheme allowed in Newton system \\
%newton\_perturbation\_merit\_factor & 0.1 & Fraction 
%of merit function value to use for perturbation in Newton system\\
%newton\_perturbation\_minimum & 1e-8 & Minimum perturbation allowed
%in 
%Newton system \\
%proximal\_perturbation & 0 & Initial perturbation to use \\
%proximal\_perturbation\_demand & no & Perturb the Newton system when 
%linear model is singular \\
%proximal\_perturbation\_first & no & Perturb the Newton system on the 
%first model only\\
%proximal\_perturbation\_shrink\_factor & 0.4 & Controls how fast the 
%perturbation shrinks\\
%proximal\_perturbation\_merit\_factor & 0.1 & Fraction of merit
%function 
%value used to decrease perturbation\\
%scale\_method & none & Scaling used for system passed to direct method\\
%\end{tabular}
%\caption{Linear Solve Options}\label{compr}
%\end{center}
%\end{table}

\section{The Nonlinear Model}\label{nonlinear}
At the nonlinear level of the algorithm, we are concerned with properties
of the algorithm affecting convergence.  These include numerical issues
related to the merit function and calculation of $H_k$ as well as
crashing and the recourse taken when a stationary point of the
merit function is encountered.  These issues are discussed in the
following subsections.  We then summarize the results and present the 
final strategies chosen.

A difficulty with the semismooth code occurs when $F$ is ill-defined because
no guarantee is made that the iterates will remain feasible with respect to 
the box $[\ell,u]$.  Such problems arise when using log functions or real
powers 
which frequently occur in applications.  Backtracking away from places 
where the function is undefined and restarting is typically sufficient for 
these models. 

\subsection{$\Phi(x^k)$ and $H_k$}
As mentioned in \secref{foundation}, the implementation of the semismooth 
algorithm for mixed complementarity problems will use the penalized 
Fischer-Burmeister merit function.
The value of $\lambda$ chosen in (\ref{CCK}) 
can have a significant impact upon the performance of the semismooth algorithm.
We note that small values of $\lambda$, say less than 0.5, should not be
used.  In the case of nonlinear complementarity problems, small values for
$\lambda$ emphasize the $\max\{0,F_i(x^k)\}\max\{0,x^k_i\}$ term of the 
penalized Fischer-Burmeister function.
This term is related to the complementarity error, but does not enforce
$F_i(x^k) \ge 0$ and $ x^k_i \ge 0$.  Since we use inexact arithmetic and 
terminate when the merit function is small, i.e. less than $ 10^{-12} $, a
small value of $\lambda$ could result in finding a point satisfying the 
termination tolerance which is not necessarily close to a solution of the 
MCP.  The default choice in the implementation is to have 
$\lambda = 0.8$ which can be changed using the {\tt chen\_lambda} 
option.

Furthermore, despite the fact that the penalized Fischer-Burmeister 
function is typically superior, there are some situations where 
the original function might be more appropriate.
Therefore, when using restarts (see Section \ref{section:restarts}) we also 
might change the merit function.  This is done by modifying the 
{\tt merit\_function} option to {\tt fischer} or {\tt chen} for the
standard and penalized Fischer-Burmeister functions.

We also note that when $\phi_1 \neq \phi_2$, there are two different
representations for $\Phi$.  This situation only occurs when both 
variables are bounded (which is the only case when $\phi_2$ is used).  In 
this case, we use 
$$
   \phi_2 (x_i - l_i, \phi_1 (u_i - x_i, - F_i (x))).
$$
An alternative is to use
$$
   -\phi_2 (u_i - x_i, \phi_1 (x_i - l_i, F_i (x))).
$$
Both of these functions are equally valid and can lead to different
sequences of iterates being generated by the algorithm.  We did not 
investigate this option further.

The calculation of $\phi(a,b)$ needs to be performed in such a way as to
minimize the effect of roundoff error.  If we have $a = 10^{-4}$ 
and $b = 10^4$
and are using a machine with 6 decimal places of accuracy, a naive
calculation of $\phi(a,b) = \sqrt{a^2 + b^2} - a - b$ would produce 
zero, meaning that we are at a solution to the problem
when in fact we are not, as the following calculation indicates:
\begin{eqnarray*}
\sqrt{10^{-8} + 10^8} - 10^{-4} - 10^4 \\
= \sqrt{10^8} - 10^{-4} - 10^4 \\
= 10^4 - 10^{-4} - 10^4 \\
= 10^4 - 10^4 \\
= 0.
\end{eqnarray*}
The actual value of $\phi(a,b)$ should be on the order of $-10^{-4}$.
A better way to calculate $\phi(a,b)$ is as follows:
\begin{enumerate}
\item If $|a| > |b|$ then $\phi(a,b) = (\sqrt{a^2 + b^2} - a) - b$.
\item Otherwise $\phi (a,b)= (\sqrt{a^2 + b^2} - b) - a$.
\end{enumerate}
This method gives a more accurate value of $\phi$.  The square root
operation needed is computed by defining $s = |a| + |b|$.  If
$s = 0$ then the value is zero, otherwise 
$s\sqrt{(\frac{a}{s})^2 + (\frac{b}{s})^2}$ is the value computed.
This eliminates overflow problems.

In order to calculate $H_k$, we use the procedure developed 
in \cite{billups:algorithms,deluca.facchinei.ea:semismooth} 
for the Fischer-Burmeister function, and a modification of the method 
in \cite{chen.chen.ea:penalized} extended for MCP models for
the penalized Fischer-Burmeister function.  
%Because of rounding
%errors, we decide that we are in the nondifferentiable area of the
%function when we are within a distance of $\epsilon$ of the region where
%$\epsilon$ is the machine precision.  The element $H_k$ we use in the
%nondifferentiable region is chosen a priori as in 
%\cite{chen.chen.ea:penalized}.  We did not experiment with choosing a
%different element in the B-subdifferential.

\subsection{Crashing}
Projected gradient crashing before starting the main algorithm can
improve the performance of the algorithm by generating a
more reasonable
starting point.  To do this, we use a technique already tested in
\cite{deluca.facchinei.ea:theoretical} and
add a new step, S.0a, to the algorithm
between S.0 and S.1.  Let $\pi_{[\ell,u]}$ be the projection of onto the
box $[\ell,u]$.  We then perform the following in this new step.

{\bf Algorithm PGC (Projected Gradient Crash)}
{\it
\begin{itemize}
\item[(S.1)] Let $j = 0$.
\item[(S.2)] Calculate $d^j = -\nabla \Psi(x^j)$.
\item[(S.3)] Let $ t_j:= \max \{ \beta^i \, | \, i = 0, 1, 2, \ldots \} $
      such that
      $$
         \Psi (\pi_{[\ell,u]}(x^j+t_j d^j)) \le \Psi (x^j) - \tau \nabla \Psi(x^j)\trt (x^j - \pi_{[\ell,u]}(x^j+t_j d^j)).
      $$
\item[(S.4)] If $t_j < \tau$ stop and set $x^0 = x^j$.
\item[(S.5)] Otherwise let $x^{j+1} = \pi_{[\ell,u]}(x^j+t_jd^j)$ and $j = j + 1$.  
Go to (S.2).
\end{itemize}
}
\noindent
In the code $\tau = 10^{-5}$ and $\beta = 0.5$;  
%These can be changed by
%setting the {\tt crash\_stepsize\_minimum} and {\tt crash\_shrink\_factor}
%respectively.  
furthermore, we only allow 10 iterations of the
projected gradient crash method.
% which can be modified with the
% {\tt crash\_iteration\_limit} option.  Crashing can be turned off with
%the option {\tt crash\_method none}.

The crashing technique presented has iterates that remain 
feasible and improve upon the initial point with respect to the merit 
function.  We believe that this is the key benefit from crashing --
all iterates remain in $[\ell,u]$.  Otherwise poor values of $x^0$ can 
frequently lead to failures in the semismooth algorithm.  In particular,
the crashing technique can significantly affect the iterates generated during
a restart.
Crashing can be turned off with the option {\tt 
crash\_method  none.}

\subsection{Stationary Points}
While stationary point termination is typically adequate for nonlinear
optimization, determining a stationary point of the residual function that
is not a zero is considered a ``failure'' by complementarity modelers.
Much theoretical work has been carried out determining the weakest
possible assumptions that can be made on the problem (and/or the
algorithm) in order to guarantee that a stationary point of the merit function
is in fact a solution of the complementarity problem.
Some of these results restrict the problem 
class considered by employing assumptions that cannot be easily 
verified for arbitrary models.  Other techniques (such as nonmonotone
linesearching) rely on a combination
of heuristics and theory, while others are entirely heuristic in nature.
The basic strategies we used to improve the reliability of 
the semismooth solver include nonmonotone linesearching and restarting.
The positive effects of these strategies have been 
demonstrated in the literature and we just present the basic idea and any
modifications made.

\subsubsection{Nonmonotone Linesearch}
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
to replace the test in step S.3 of Algorithm BSM with the nonmonotone
test:
$$
         \Psi (x^k+t_k d^k) \le R^k + t_k \nabla \Psi(x^k)\trt d^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 reference value choice.  We begin by letting
$\{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. The value of $\kappa$ defaults to 1 in the
code and can be modified with the {\tt nms\_initial\_reference\_factor} 
option;
% and is used to determine the initial set of acceptable 
%merit function values.  
$\kappa=1$ indicates that we are not going to 
allow the merit function to increase beyond its initial value.  
%The code uses 
%$m = 4$, which can be changed with the {\tt nms\_memory\_size} option.

Having defined the values of $\{M_1,\ldots,M_m\}$ (where the code by default
uses $m=4$), we can now calculate a reference
value.  We must be careful when we allow gradient steps in the code.
Assuming that $d^k$ is the Newton direction (or a least squares solution
to the Newton system in the presence of singularity, see \ref{least}),
we define $i_0 = \mbox{argmax }M_i$ and $R^k = M_{i_0}$.  
After the nonmonotone 
linesearch rule above finds $t_k$, we update the memory so
that $M_{i_0} = \Psi(x^k+t_kd^k)$. That is, we remove an element from 
the memory having the largest merit function value.

When we decide to use a gradient step, 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 to say that 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 $m$ iterations,
where $m$ is set by the {\tt nms\_mstep\_frequency} option.  Currently,
we use this strategy only used in a restart.

\subsubsection{Restarting}\label{section:restarts}
The rules for nonmonotone linesearching and crashing are extremely 
useful in practice, but do not preclude convergence to a nonoptimal
stationary point.  One observation relevant for complementarity solvers is 
that we know a priori the optimal value of the merit function at a solution 
if one exists. If the 
code detects that the current iterate is a stationary point that is not a 
solution, a recovery strategy can be invoked.  One successful technique is the
restart strategy introduced in PATH where the recovery mechanism involves 
starting over from the user supplied starting point with a different set of 
options, thus leading to a different sequence of iterates being investigated.
For the semismooth algorithm we use the restarts defined in 
Table \ref{restarts}.
In addition, the first restart of the semismooth code will also include 
the option {\tt chen\_lambda 0.95} if no crash iterations were performed
in its first run; otherwise, we would waste computational
resources by generating the exact same sequences of iterates as the
first attempt.

\begin{table}
\begin{center}
\caption{Restart Definitions}\label{restarts}
\begin{tabular}{ll}
Restart Number & Parameter Values \\
\hline
1 & crash\_method none \\
\hline
2 & chen\_lambda 0.95 \\
  & crash\_method none \\
  & nms\_mstep\_frequency 4 \\
  & nms\_initial\_reference\_factor 5 \\
\hline
3 & crash\_method none \\
  & merit\_function fischer \\
  & nms\_mstep\_frequency 1 \\
  & nms\_initial\_reference\_factor 1 \\
\end{tabular}
\end{center}
\end{table}

Restarts should be applicable to general models and not written to overcome 
difficulties with a specific model.  That is, if we were to use the restart 
definition as the default, we should still solve most problems in the test 
set.  We present in
Table \ref{restart} the numbers of failures on a subset of GAMSLIB and
MCPLIB models when the 
particular restart options were used.  Note that the restarts (0, 1 and 2) 
using the penalized Fischer-Burmeister function given in equation (\ref{CCK}) 
outperform the one using the standard
Fischer-Burmeister function defined in equation (\ref{FB}).

\begin{table}
\begin{center}
\caption{Restart Performance on GAMSLIB and MCPLIB Problems}\label{restart}
\begin{tabular}{lcc}
               & \multicolumn{2}{c}{Failures} \\
Restart Number & GAMSLIB & MCPLIB \\
\hline
0 (first run) & 3 & 63 \\
1 & 3 & 63 \\
2 & 3 & 64 \\
3 & 12 & 76
\end{tabular}
\end{center}
\end{table}

\section{Computational Results}
The LUSOL \cite{gill.murray.ea:maintaining} 
sparse factorization routines contained in the 
MINOS \cite{murtagh.saunders:minos} nonlinear programming solver were used
for factorization purposes.
All of the linear algebra and other basic mechanisms are exactly the same
for both SEMI and PATH 4.x.  Therefore, the comparison made 
is as close 
to a true comparison of the algorithms as we can make.  We note that 
PATH 4.x is much more mature than the semismooth implementation.  Both 
codes are continually being improved when deficiencies are uncovered.

Finally, we remark that the factorization routine is a key component of
the implementation.  For the semismooth algorithm, rank-1 updates to the 
matrix are not required.
Therefore, factorization routines other than those
provided by LUSOL might be more effective for SEMI.  
However, we sacrifice some potential speed gains in preference for having
a consistent choice for the direct method among the codes tested.

To test the standard algorithm, we ran the code on all of the problems in
the GAMSLIB and MCPLIB \cite{dirkse.ferris:mcplib} suites of test problems.
The following tables report the number of successes, followed by the number 
of failures in parenthesis.  In order to test the reliability of the 
iterative method, we ran all of the models using only the iterative LSQR 
technique to calculate the Newton direction.  For comparison purposes, we 
also show the performance of PATH 4.x from \chapref{path}
on the same problems. These latter two results are given in the columns
labeled ``Iterative'' and ``PATH 4.x'' respectively.  We split the results 
into those that allow restarts and those that do not.  

These results indicate that while SEMI is not
as reliable as PATH 4.x, it does well and is robust on the GAMSLIB and
MCPLIB problem sets.  We remark that the \texttt{box} model accounts for 
over 80 failures in the MCPLIB test set.  This particular test problem 
modifies a parameter within a loop and solves each successive problem using 
the solution from the previous iteration as a starting point.  Therefore, a 
failure to solve one particular problem leads to extreme difficulty solving 
the next.  Hence, we have a large number of accumulated failures for this
problem.  We 
could argue for the removal of the \texttt{box} model from the test set.
However, for completeness we report results for the \texttt{box} model.

The results also indicate that the iterative method is robust.
Since we are only using iterative techniques in this code, the memory
requirements are much less.  42 of the additional failures in the 
``Iterative'' column occur in the {\tt asean9a}, {\tt denmark}, and 
{\tt opt\_cont511} models in which the time limit was encountered.  We also 
note that the restart heuristic significantly improves the robustness of both 
SEMI and PATH 4.x on MCPLIB.

\begin{comment}
The results on the NETLIB test set are disappointing, but expected.  The 
BD-regularity assumption made in Theorem \ref{Convergence} does not typically 
hold for the optimality conditions of linear programs, as the Jacobian of
$F$ is only a $P_0$-matrix.  
%Note that if we have a linear complementarity 
%problem defined by a P-matrix, then Proposition \label{CharP} can be used
%to show that the solution, which is exists and is unique, is a BD-regular 
%solution to $\Phi$. 
Algorithms have been proposed for standard linear complementarity problems 
defined by a positive semidefinite matrix \cite{fischer:local}, however, we 
did not investigate these methods.
\end{comment}

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

\begin{table}
\centering
\caption{Comparison of SEMI and PATH 4.x on MCPLIB}
\begin{tabular}{l|ccc|ccc}
Problem & \multicolumn{3}{c|}{With Restarts} &
          \multicolumn{3}{c}{Without Restarts} \\
        & SEMI & Iterative & PATH 4.x & SEMI & Iterative & PATH 4.x \\
\hline
asean9a &1 (0) &	0 (1) &	1 (0) &	1 (0) &	0 (1) &	1 (0) \\
badfree & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
bai\_haung & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
bert\_oc & 4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) \\
bertsekas & 6 (0) &	6 (0) &	6 (0) &	6 (0) &	6 (0) &	6 (0) \\
billups & 0 (3) &	0 (3) &	0 (3) &	0 (3) &	0 (3) &	0 (3) \\
bishop & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
box & 277 (84) &	280 (81) &	361 (0) &	270 (91) &	271 (90) &	361 (0) \\
bratu & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
cammcf & 0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) \\
carbon & 6 (1) &	6 (1) &	6 (1) &	6 (1) &	6 (1) &	6 (1) \\
cgereg & 17 (2) &	18 (1) &	19 (0) &	16 (3) &	18 (1) &	19 (0) \\
choi & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
colvdual & 4 (0) &	4 (0) &	4 (0) &	3 (1) &	3 (1) &	2 (2) \\
colvnlp & 6 (0) &	6 (0) &	6 (0) &	6 (0) &	6 (0) &	6 (0) \\
congest & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
cycle & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
degen & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
denmark & 38 (0) &	0 (40) &	38 (0) &	28 (10) &	0 (40) &	38 (0) \\
dirkse1 & 0 (2) &	0 (2) &	1 (1) &	0 (2) &	0 (2) &	1 (1) \\
duopoly & 0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) \\
eckstein & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
ehl\_k40 & 3 (0) &	3 (0) &	3 (0) &	3 (0) &	3 (0) &	2 (1) \\
ehl\_k60 & 3 (0) &	3 (0) &	3 (0) &	3 (0) &	3 (0) &	2 (1) \\
ehl\_k80 & 3 (0) &	3 (0) &	3 (0) &	3 (0) &	3 (0) &	3 (0) \\
ehl\_kost & 3 (0) &	3 (0) &	3 (0) &	3 (0) &	3 (0) &	3 (0) \\
electric & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	0 (1) \\
endog & 21 (2) &	0 (2) &	23 (0) &	20 (3) &	0 (2) &	23 (0) \\
eppa & 8 (0) &	8 (0) &	8 (0) &	8 (0) &	8 (0) &	8 (0) \\
eta2100 & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
exemptions & 182 (4) &	182 (4) &	186 (0) &	182 (4) &	182 (4) &	186 (0) \\
explcp & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
exros & 5 (0) &	5 (0) &	5 (0) &	5 (0) &	5 (0) &	2 (3)
\end{tabular}
\end{table}

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

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

\begin{table}
\centering
\caption{Comparison of SEMI and PATH 4.x on MCPLIB (cont.)}
\begin{tabular}{l|ccc|ccc}
Problem & \multicolumn{3}{c|}{With Restarts} &
          \multicolumn{3}{c}{Without Restarts} \\
        & SEMI & Iterative & PATH 4.x & SEMI & Iterative & PATH 4.x \\
\hline
ralph & 7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) \\
renger & 0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) \\
romer & 2 (0) &	2 (0) &	2 (0) &	2 (0) &	2 (0) &	1 (1) \\
romer2 & 0 (7) &	0 (7) &	1 (6) &	0 (7) &	0 (7) &	0 (7) \\
romer3 & 1 (0) &	0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) \\
runge & 7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) \\
scarfanum & 4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) \\
scarfasum & 4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) \\
scarfbnum & 2 (0) &	2 (0) &	2 (0) &	2 (0) &	2 (0) &	2 (0) \\
scarfbsum & 2 (0) &	2 (0) &	2 (0) &	2 (0) &	2 (0) &	2 (0) \\
shubik & 46 (1) &	46 (1) &	47 (0) &	38 (9) &	38 (9) &	43 (4) \\
simple-ex & 1 (0) &	1 (0) &	0 (1) &	0 (1) &	0 (1) &	0 (1) \\
simple-red & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
spill1 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
spill2 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	0 (1) \\
spill3 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	0 (1) \\
spill4 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
spill5 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
sppe & 3 (0) &	3 (0) &	3 (0) &	3 (0) &	3 (0) &	3 (0) \\
tinloi & 64 (0) &	64 (0) &	64 (0) &	64 (0) &	64 (0) &	56 (8) \\
tinloilp & 0 (4) &	0 (4) &	4 (0) &	0 (4) &	0 (4) &	4 (0) \\
tinsmall & 64 (0) &	64 (0) &	64 (0) &	64 (0) &	64 (0) &	63 (1) \\
titan & 1 (1) &	1 (1) &	1 (1) &	1 (1) &	1 (1) &	1 (1) \\
tobin & 4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) \\
tqbilat & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
trade12 & 2 (0) &	2 (0) &	2 (0) &	2 (0) &	2 (0) &	2 (0) \\
trafelas & 1 (1) &	0 (2) &	2 (0) &	1 (1) &	0 (2) &	2 (0) \\
trig & 2 (1) &	2 (1) &	3 (0) &	0 (3) &	0 (3) &	3 (0) \\
uruguay & 4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) &	4 (0) \\
venables & 2 (0) &	2 (0) &	2 (0) &	2 (0) &	2 (0) &	2 (0) \\
vonthmcf & 0 (1) &	1 (0) &	1 (0) &	0 (1) &	1 (0) &	1 (0) \\
water & 3 (0) &	3 (0) &	3 (0) &	3 (0) &	3 (0) &	2 (1) \\
xu1 & 7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) \\
xu2 & 7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) \\
xu3 & 7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) \\
xu4 & 7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) \\
xu5 & 7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) &	7 (0) \\
\hline
Total & 1081 (151) &	1024 (189) &	1207 (24) &	1034 (198) &	987 (226) &	1171 (60)
\end{tabular}
\end{table}

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

\begin{table}
\centering
\caption{Comparison of SEMI and PATH 4.x on NETLIB (cont.)}
\begin{tabular}{l|ccc|ccc}
Problem & \multicolumn{3}{c|}{With Restarts} &
          \multicolumn{3}{c}{Without Restarts} \\
        & SEMI & Iterative & PATH 4.x & SEMI & Iterative & PATH 4.x \\
\hline
israel & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
kb2 & 1 (0) &	1 (0) &	1 (0) &	0 (1) &	1 (0) &	1 (0) \\
lotfi & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
nesm & 0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) \\
perold & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
pilot4 & 0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) &	0 (1) \\
sc105 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
sc205 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
sc50a & 0 (1) &	1 (0) &	1 (0) &	0 (1) &	1 (0) &	1 (0) \\
sc50b & 0 (1) &	1 (0) &	1 (0) &	0 (1) &	1 (0) &	1 (0) \\
scagr25 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
scagr7 & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
scfxm1 & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	0 (1) &	1 (0) \\
scfxm2 & 1 (0) &	0 (1) &	1 (0) &	1 (0) &	0 (1) &	1 (0) \\
scfxm3 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
scorpion & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
scrs8 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
scsd1 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
scsd6 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
scsd8 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
sctap1 & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
sctap2 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
sctap3 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
share1b & 1 (0) &	1 (0) &	1 (0) &	0 (1) &	1 (0) &	1 (0) \\
share2b & 1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) &	1 (0) \\
shell & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
ship04l & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
ship04s & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
ship08l & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
ship08s & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
ship12l & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
ship12s & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
stair & 1 (0) &	1 (0) &	1 (0) &	0 (1) &	1 (0) &	1 (0) \\
stocfor1 & 1 (0) &	1 (0) &	1 (0) &	0 (1) &	1 (0) &	1 (0) \\
stocfor2 & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
tuff & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
wood1p & 0 (1) &	0 (1) &	1 (0) &	0 (1) &	0 (1) &	1 (0) \\
\hline
Total & 22 (53) &	22 (53) &	66 (9) &	14 (61) &	17 (58) &	66 (9)
\end{tabular}
\end{table}
\end{comment}

We currently do not have any results on very large problems, but believe that 
based on the evidence, the code will scale well to the larger problems.  In
particular, the iterative version of the SEMI requires
significantly less memory than PATH 4.x, allowing the possibility of
solving huge models. 
The current drawback of the semismooth code is the time taken by LSQR
to solve the linear systems.  We designed the code for robustness, and
therefore chose parameters in the code to enhance reliability.  
This choice results in PATH 4.x being much faster than the semismooth 
code.  However, we believe that many of the standard techniques for 
improving iterative linear equation solvers are applicable in the context 
of the semismooth algorithm.

\section{Summary}
We have shown that a code based on semismooth algorithm can be implemented 
as a robust MCP solver.  Particular care needs to be taken in implementing
the evaluation of $\Phi$ and $\Psi$, and in treating numerical issues
related to constructing the Newton step.  The resulting SEMI code is
robust and has great potential in the very large scale setting as indicated 
by the numerical results.

We remark that SEMI is the first implementation of a semismooth algorithm 
for solving mixed complementarity problems that is available in the AMPL 
and GAMS modeling languages and the MATLAB and NEOS tools.  We also note that
an earlier version of this chapter appeared in 
\cite{munson.facchinei.ea:semismooth}.

