
\chapter{Preprocessing}\label{chap:preprocessing}

% This chapter is an update to sections of the following papers:
%   ``Preprocessing Complementarity Problems''
%   ``Case Studies in Complementarity: Improving Model Formulation''
%       - Diagnostic information
% Added the documentation on the options effecting the code.

The benefits of preprocessing have long been known to the 
linear \cite{andersen.andersen:presolving,brearley.mitra.ea:analysis} and 
mixed integer \cite{savelsbergh:preprocessing} programming communities, yet 
have not previously been studied from a complementarity perspective.  
In many cases, the preprocessor simply uncovers structure from a model that 
may not have been explicitly communicated.  The reasons for this are varied.
The modeler may be completely unaware of the structure, or unwilling or 
unable to pass the known structure to a solver.  In other cases, the 
structure arises from simplifications due to particular data values.  For 
example, if a general inequality constraint can be reduced to use only one 
variable, the solver might prefer to treat it as a bound constraint.

In general, the purpose of a preprocessor is to reduce the size and 
complexity of a model to achieve improved performance by the main algorithm.
For example, by fixing variables, some nonlinear constraints may become 
linear, and smaller linear systems are typically solved with less memory and 
in less time.  Another benefit of the analysis performed is the detection of 
model infeasibilities.  Most codes continue attempting to solve a model until 
a time or iteration limit is reached.  In cases where the model has no 
solution, this is obviously wasteful.  A preprocessor can sometimes process 
the model enough to determine (and report) the source of the difficulty.

The preprocessor for complementarity problems developed in this chapter uses 
the box constrained variational inequality formulation 
\cite{karamardian:complementarity,robinson:generalized*1} of the mixed 
complementarity problem:
\begin{eqnarray}
0 \in F(x) + N_{[\ell,u]}(x), \label{eqn:model}
\end{eqnarray}
where $F: \real^n \rightarrow \real^n$ is continuously differentiable and 
\begin{eqnarray*}
\ell \in \{ \real\cup\{-\infty\}\}^n\\
u \in \{ \real\cup\{+\infty\}\}^n
\end{eqnarray*}
with $[\ell,u]$, the Cartesian product of the closed, not necessarily compact, 
intervals $[\ell_i,u_i]$.  Recall that the normal cone 
\cite{rockafellar:convex} to the box $[\ell,u]$ can be determined as a 
Cartesian product, namely
\[
N_{[\ell,u]}(x) = \prod_{i=1}^{n} N_{[\ell_i,u_i]}(x_i).
\]
Each component in this product depends on $x_i$, $\ell_i$ and
$u_i$ but is either $\real$, $\real_+$, $\real_-$, $\{ 0 \}$, or 
$\emptyset$.  In particular,
\[
N_{[\ell_i,u_i]}(x_i) = 
\left\{
  \begin{array}{ll} 
    \real & \text{if } \ell_i = x_i = u_i \\
    \real_- & \text{if } \ell_i = x_i < u_i\\
    \{0\} & \text{if } \ell_i < x_i < u_i \\
    \real_+ & \text{if } \ell_i < x_i = u_i \\
    \emptyset & \text{otherwise.}
  \end{array}
\right.
\]
See \defref{normal} for a specification of the normal cone for a general
closed convex set.

% As an example, suppose
% \[ [\ell,u] = [-2,\infty) \times \{ 3 \} \subseteq \real^2 . \]
% In this case, 
% \[ N_{[\ell,u]}(x) = 
% \left\{
% \begin{array}{cl} 
% \real_- \times \real & \mbox{if } x_1 = -2, x_2 = 3 \\
% \{ 0 \} \times \real & \mbox{if } x_1 > -2, x_2 = 3 \\
% \emptyset & \mbox{otherwise,}
% \end{array}\right.
% \]
% and hence \eqnref{model} is the mixed system of inequalities and
% equations
% \begin{eqnarray*}
% F_1(x) & \geq & 0 \\
% x_1 & \geq & -2 \\
% (x_1 + 2) F_1(x) & = & 0 \\
% x_2 & = & 3
% \end{eqnarray*}
% Note that $F_2$ is irrelevant, and so the example reduces to the 
% standard complementarity problem 
% \[ 
% F_1(z-2,3) \geq 0 \perp z \geq 0
% \] 
% in the single variable $z$.

The preprocessor works upon two equivalent representations of the same 
complementarity problem.  To understand the basic methodology developed, 
consider a standard convex quadratic programming problem:
\begin{eqnarray}
\begin{array}{lc}
\mbox{min} & \frac{1}{2}x^TQx + c^Tx \\
\mbox{s.t.} & Ax \ge b \\
& x \ge 0,
\end{array}\label{eqn:qp}
\end{eqnarray}
where $Q \in \real^{n \times n}$ is a symmetric positive semidefinite matrix, 
$A \in \real^{m \times n}$, $b \in \real^m$, and $c \in \real^n$. The
quadratic program in \eqnref{qp}
can be posed as a variational inequality in one of two ways.  First, when dual
variables, $\lambda$, are introduced, the complementary slackness conditions 
for quadratic programs form the box constrained variational inequality:
\begin{eqnarray}
0 \in 
\left[\begin{array}{cc} Q & -A^T \\ A & 0 \end{array}\right] 
\left[\begin{array}{c} x \\ \lambda \end{array}\right] + 
\left[\begin{array}{c} c \\ -b \end{array}\right] +
\left[\begin{array}{c} N_{\real^n_+}(x) \\ N_{\real^m_+}(\lambda)\end{array}\right].
\label{eqn:qpcomp}
\end{eqnarray}
Alternatively, the first order conditions can be succinctly written as
a polyhedrally constrained variational inequality:
\begin{eqnarray}
0 \in Qx + c + N_C(x), \label{eqn:qpfirst}
\end{eqnarray}
where $C = \{x \mid Ax \ge b, x \ge 0\}$.  Since $C$ is a geometric object, 
a computationally attractive algebraic representation can be chosen for $C$. 
Exploiting two representations of the mixed complementarity problem
(analogous to \eqnref{qpcomp} and \eqnref{qpfirst} for quadratic programs) 
is a key concept in the preprocessor.

The complementarity problem is communicated as the box constrained 
variational inequality \eqnref{model}.  Therefore, any polyhedral 
structure will need to be recovered from the problem description.  
Once recovered, the preprocessor automatically rewrites the problem 
in the form
\begin{equation}
0 \in H(x) + N_C(x) \label{eqn:model-red}
\end{equation}
for a polyhedral set $C$.  
Both \eqnref{model} and \eqnref{model-red} are used in distinct phases 
of the preprocessor.  
For example, the general inequalities in the set $C$ can be used 
to modify the bounds $\ell_i$ and $u_i$ on a variable $x_i$.  In particular,
if $x_i = \ell_i = u_i$ then $N_{[\ell_i,u_i]}(x_i) = \real$.  Hence, 
fixing a variable $x_i$ means that the corresponding constraint
\[
0 \in F_i(x) + N_{[\ell_i,u_i]}(x_i) \iff 0 \in F_i(x) + \real
\]
is trivially satisfiable.  Thus, preprocessing in the complementarity 
case attempts to fix variables and thus remove constraints.

\secref{polyhedral} begins by detailing the process used
to uncover and exploit polyhedral structure in an MCP.  The general idea 
is to reformulate \eqnref{model} in a form similar to \eqnref{qpfirst},
with a general polyhedral set $C$ replacing $[\ell,u]$.  The representation 
of the set $C$ can then be modified by either removing constraints or 
bounding variables.  When converted back to a mixed complementarity problem 
a reduction in the number of variables is realized.  Note that the process 
developed in this section recovers most checks done by traditional linear 
programming codes \cite{andersen.andersen:presolving} when given the 
complementary slackness necessary and sufficient conditions for linear 
programs, but is applicable to a larger class of problems.

Further reductions to the MCP can be made by utilizing information about 
$F$ and its Jacobian, $F'$, as developed in \secref{function}.  In 
particular, the range of $F$ is used to eliminate variables from the model, 
row and column duplicates can also be removed, and by detecting special block 
structure, a sequence of smaller problems can be solved to find an answer 
to the original problem.

Both phases are incorporated into a complete preprocessor for 
mixed complementarity problems in \secref{preproc-results}.  
Computational results for some test problems are presented indicating 
the success of the procedure outlined.

More information about the problem must be provided to the preprocessor 
than is necessary to solve it.  The basic requirement is a listing of the 
linear and nonlinear elements in the Jacobian of $F$, which are conveyed
to the solver using the \texttt{Presolve\_Interface} from 
\chapref{interfaces}.  
This knowledge is sufficient to find and utilize special structures.  The AMPL 
\cite{fourer.gay.ea:ampl} and GAMS \cite{brooke.kendrick.ea:gams} 
environments, and the \texttt{pathlcp} and \texttt{semilcp} functions in 
MATLAB already provide this information.  Users of other interfaces, such as 
NEOS \cite{ferris.mesnier.ea:neos} and the \texttt{pathmcp} and 
\texttt{semimcp} functions in MATLAB, will need to develop the appropriate 
routines.
Some checks in \secref{function} based on the nonlinear functions 
need to know the range of $F$ over $[\ell,u]$.  Routines to calculate these 
intervals are not currently provided by any of the interfaces.

\section{Polyhedral Constraints}\label{sec:polyhedral}
The first stage of the preprocessor detects polyhedral structure in a mixed 
complementarity problem and exploits it by transforming the source problem 
into a model of lower dimension where $C$ is the intersection 
of a closed product of intervals and a polyhedral set.  The representation 
of $C$ is then modified by removing constraints and changing bounds with
the resultant MCP having fewer variables.  After the preprocessed model has 
been solved, a solution to the original MCP is recovered with a postsolve step.

\subsection{Presolve}
Polyhedral structure can be exploited when given a special type of 
complementarity problem.  Suppose the variables can be split into $(x,y)$ and
\eqnref{model} has the form:
\begin{eqnarray} 
0 \in \left[\begin{array}{c} F(x) - A^Ty \\ Ax - b\end{array}\right] +
\left[\begin{array}{c} N_X(x) \\ N_Y(y)\end{array}\right],\label{eqn:bvi}
\end{eqnarray}
where $F: \real^n \rightarrow \real^n$ is continuously differentiable,
$A \in \real^{m \times n}$, $b \in \real^m$, $X \subseteq \real^n$ is a 
Cartesian product of closed intervals, and $Y \subseteq \real^m$ is a 
Cartesian product of $\real$, $\real_+$, or $\real_-$.  Note that if $Y_i = \{0\}$ 
then $y_i$ can be fixed at zero and the corresponding 
$A_{i,\cdot}x - b_i$ removed.  Further, if $Y_i = [\ell_i, \infty)$ or 
$Y_i = (-\infty, u_i]$ for some finite $\ell_i$ or $u_i$ then an appropriate 
change of variables, possibly adding constant vectors to $F(x)$ and $b$,
replaces $Y_i$ with $\real_+$ and $\real_-$ respectively.

\begin{theorem}\label{thm:robinson}(Propositions 1 and 2 of 
\cite{robinson:reduction})
Under the assumptions placed on $X$, $Y$, and the structure of the problem 
given above the following hold:
\begin{itemize}
\item[(a)] If $(\bar{x}, \bar{y})$ solves \eqnref{bvi} then $\bar{x}$ is a solution
to
\begin{eqnarray}
0 \in F(x) + N_{X \cap \{x \mid b - Ax \in Y^o\}}(x),\label{eqn:pvi}
\end{eqnarray}
where $Y^o$ denotes the polar cone of $Y$ which is defined as
\[
Y^o := \{ y \mid \ip{y}{\bar{y}} \le 0, \;\forall\bar{y} \in Y\}.
\]
\item[(b)] If $\bar{x}$ solves \eqnref{pvi} then the linear optimization problem
\begin{eqnarray}
\begin{array}{cc} \min_{y \in Y} & \ip{A\bar{x} - b}{y} \\
\mbox{s.t.} & 0 \in F(\bar{x}) - A^Ty + N_X(\bar{x}).
\end{array}\label{eqn:opt}
\end{eqnarray}
has a nonempty solution set.  Further, for any $\bar{y}$ 
solving \eqnref{opt}, $(\bar{x}, \bar{y})$ solves \eqnref{bvi}.
\end{itemize}
\end{theorem}
%{\bf Proof.} The first assertion follows directly from Propositions 1 and
% 2 of \cite{robinson:reduction}.  The first part of the second assertion 
% follows from the comments following Proposition 1 in 
% \cite{robinson:reduction}, while the second applies Proposition 1.

\thmref{robinson} provides the machinery used by the first stage of the
preprocessor.  The Jacobian matrix, $F'$, is stored in both row- and 
column-oriented data structures.  Utilizing the information provided about the 
types of the elements in the Jacobian, a row and column possessing the 
necessary skew symmetric structure of \eqnref{bvi} can be quickly identified.
\thmref{robinson} is then applied to this single row and column to
create a problem of the form \eqnref{pvi}.  The polyhedral set, 
$X \cap \{x \mid b - Ax \in Y^o\}$,  is then checked
for possible reductions, that is whether the general constraint,
$b - Ax \in Y^o$, can be moved into the bound constraint, $X$.  The
new set $\tilde{X} \cap \{x \mid b - Ax \in \tilde{Y}^o\}$ is identical to
$X \cap \{x \mid b - Ax \in Y^o\}$ but the MCP recovered using 
\thmref{robinson} on the reduced model is typically simpler.
The identification and modification continues until no further 
simplifications are made.
Note that \thmref{robinson} is only applied to a small number of constraints 
at a time during preprocessing while \cite{robinson:reduction} uses the
machinery to ready a problem for solution by a polyhedrally constrained
variational inequality solver \cite{sellami.robinson:implementation}.  
Automatically finding a set of polyhedral constraints with 
maximum size from \eqnref{model} is a harder problem and is not 
considered here.

Information about any modifications performed are pushed onto a stack.  The 
stack is a convenient data structure with two basic operations: pushing an 
element onto the top and popping an element from the top.  Changes are pushed 
onto the stack in the order performed and are popped off the stack in the 
reverse order during the postsolve.

The complete algorithm for the first phase of the preprocessor is as follows:

{\bf Algorithm Presolve}
{\it
\begin{itemize}
\item[(S.1)] Mark all rows and columns with the skew symmetric structure
as {\em eligible} candidates, excluding any rows complementary to a variable with
finite lower and upper bounds.
\item[(S.2)] Using some ordering, pick one of the candidates and transform the 
problem into a polyhedrally constrained equation using \thmref{robinson}.
\item[(S.3)] Analyze the polyhedral set and modify the representation as
detailed below.  Push any changes on top of the stack.
\item[(S.4)] Transform the modified problem back to box constrained form.
\item[(S.5)] Repeat steps (S.2)--(S.4) until there are no reductions possible.
\end{itemize}
}

\noindent
The implementation performs all simple reductions (\secref{simple} and
\secref{doubleton}) first.
Once all of these are completed, forcing constraints (\secref{force}) and 
redundant rows (\secref{redundant}) are checked.  In a nonlinear model, 
additional rows and columns can become linear when variables are fixed.  
Therefore, after all tests are completed on the current list of eligible 
candidates, another pass is made through the Jacobian to mark new eligible 
rows and columns which are checked using (S.2)--(S.5) of Algorithm Presolve.  
When no new eligible rows and columns are created the process stops.

\subsubsection{Simple Reductions}\label{sec:simple}
The simplest reduction made is when an eligible row contains zero elements.  
This situation corresponds to the case where the polyhedral set in the 
transformed problem is:
\begin{eqnarray}
        X \cap \{x \mid b \in Y^o\}\label{set:zero}.
\end{eqnarray}
If \setref{zero} is empty, then the original problem has no solution.  
Otherwise, the polyhedral component is irrelevant, because it does not
impose any constraints on $x$, and \setref{zero} can be replaced with:
\[
        X \cap \{x \mid b \in \{0\}^o\}.
\]
Note that $\{0\}^o = \real$, ensuring the constraint $b \in \{0\}^o$
is always satisfied.  When 
transformed back to the original space, the corresponding 
multiplier is fixed at 0 and 
removed from the problem, resulting in a reduction of one variable and the
corresponding constraint.

Another simple reduction occurs when the eligible row contains a single
element.  The polyhedral set in this case is:
\begin{eqnarray}
        X \cap \{x \mid b - ax_i \in Y^o\}\label{set:one}.
\end{eqnarray}
Since $Y$ is $\real_+$, $\real_-$, or $\real$, the constraint will be either
$b - ax_i \le 0$, $b - ax_i \ge 0$, or $b - ax_i = 0$ respectively.
Each of these imply simple bounds on $x_i$, which can be explicitly 
incorporated in $X$.  Therefore, \setref{one} is replaced by
\[
        \tilde{X} \cap \{x \mid b - ax_i \in \{0\}^o\},
\]
where $\tilde{X}$ includes the tightened bounds on $x_i$.  This modification
results in a reduction of at least one variable.

\subsubsection{Doubleton Rows}\label{sec:doubleton}
Doubleton rows having the skew symmetric property can also be preprocessed.
Consider a row, $i$, which is an equation of the form
\[
    ax_j + bx_k = c
\]
where either $x_j$ or $x_k$ is a column singleton.  Assume $x_k$ is the
column singleton.  If row $k$ also has the skew symmetric property, 
then $X$ can be modified by changing the bounds on $x_j$ to 
make $x_k$ free.  Immediately following this change, row $k$ (which
must be a singleton row complementary to the free variable $x_k$) is
preprocessed out of the model using the singleton check described above.

\subsubsection{Forcing Constraints}\label{sec:force}
Forcing constraints are constraints for which, given the bounds on the 
variables, there is exactly one feasible point.  Once it is known that
only one solution is possible, all variables appearing in the constraint can 
be fixed, potentially leading to a large reduction in problem size.

Let the polyhedral constraint be written in the form:
\begin{eqnarray}
        X \cap \{x \mid b - a^Tx \in Y^o\}\label{set:force}.
\end{eqnarray}
Without loss of generality, assume that $Y = \real_+$ which means
that $Y^o = \real_-$.  Then \setref{force} can then be explicitly stated as:
\[
        X \cap \{x \mid b \le a^Tx\}.
\]
Using $X$, bounds, $\underline{a}$ and $\bar{a}$, can be implied such that
$\underline{a} \le a^Tx \le \bar{a}$ for all $x \in X$.
The ranges are determined as follows:
\begin{eqnarray*}
\underline{a} = \sum_{\{i \mid a_i > 0\}} a_i\ell_i + \sum_{\{i \mid a_i < 0\}} a_iu_i \\
\bar{a} = \sum_{\{i \mid a_i > 0\}} a_iu_i + \sum_{\{i \mid a_i < 0\}} a_i\ell_i, \\
\end{eqnarray*}
where $\ell_i$ and $u_i$ are the lower and upper bounds on variable $i$ 
respectively.  If $\bar{a} < b$ the problem is infeasible.  Otherwise, if 
$\bar{a} = b$, there is only one feasible point and the values of 
the variables are fixed at $u_i$ for all $i$ with $a_i > 0$ and $\ell_i$ 
for all $i$ with $a_i < 0$.  Set \setref{force} is then replaced with:
\begin{eqnarray}
        \tilde{X} \cap \{x \mid b - a^Tx \in \{0\}^o\},
\end{eqnarray}
where $\tilde{X}$ contains the fixed variables.  The net result is that
the forcing constraint and a number of variables are removed from the 
original problem.

\subsubsection{Redundant Rows}\label{sec:redundant}
Redundancy in the Jacobian matrix can cause difficulty for many algorithms.
Therefore, it is advantageous to remove as much redundancy as possible.
The algorithm given in \cite{tomlin.welch:finding} is used to identify 
duplicate rows.  All eligible constraints are checked simultaneously with the 
algorithm.  After finding two duplicate rows, any inconsistencies are
uncovered (meaning that the model is unsolvable) and one of the constraints is
removed wherever possible.
Without loss of generality, let the constraint set be written as
\begin{eqnarray}
        X \cap \left\{{x \left|\begin{array}{c} b - ax \in Y_1^o \\
 c - ax \in Y_2^o\end{array}\right.}\right\}\label{set:redundant}.
\end{eqnarray}
Several cases are presented in \tabref{redundant} along with the associated 
actions taken.  The other cases are symmetric to those given in the table.

\begin{table}
\centering
\caption{Redundant Rows Cases}\label{tab:redundant}
\begin{tabular}{lp{10.5cm}}
Case & Action \\
\hline
$Y_1 = \real$ and $Y_2 = \real$   & If $b = c$ remove one of the constraints.
                                Otherwise the problem is infeasible.\\
$Y_1 = \real$ and $Y_2 = \real_+$ & If $b \ge c$ remove the inequality
constraint.
                               Otherwise the problem is infeasible.\\
$Y_1 = \real_+$ and $Y_2 = \real_+$ & If $b \ge c$ remove the constraint associated 
                                with $Y_2$. 
				Otherwise remove the $Y_1$ 
                                constraint.\\
$Y_1 = \real_-$ and $Y_2 = \real_+$ & If $b < c$, the problem is infeasible. 
                                Otherwise if $b = c$ make one of the 
                                constraints an equation and remove the other.
                                Otherwise, it is a range constraint; 
                                nothing is done by the preprocessor.\\
\end{tabular}
\end{table}

\subsubsection{Extensions}
The requirements in \thmref{robinson} can be slightly weakened.
Let $D \in \real^{n \times n}$ be a positive diagonal matrix.  Then the
following form will suffice instead of \eqnref{bvi}:
\begin{eqnarray}
0 \in \left[\begin{array}{c}F(x) - DA^Ty\\ Ax - b\end{array}\right] + 
\left[\begin{array}{c}N_X(x)\\N_Y(y)\end{array}\right]\label{eqn:bvi2}
\end{eqnarray}
because \eqnref{bvi2} can be reduced to \eqnref{bvi} by applying a
diagonal row scaling of 
$\left[\begin{array}{cc}D & 0 \\ 0 & 1\end{array}\right]^{-1}$
and recalling that the normal cone does not change under multiplication by a 
positive diagonal matrix.

Furthermore, free variables imply that the corresponding function is an
equation because $N_\real(\cdot) \equiv \{0\}$.  Therefore, the equation can
be negated to obtain the required structure.

\subsection{Postsolve}
Once the algorithm has solved the preprocessed model, all of the presolve 
steps must be undone in the reverse order to recover a solution to the 
original model.  The stack of presolve records is used for this purpose.  
The following algorithm is performed:

{\bf Algorithm Postsolve}
{\it
\begin{itemize}
\item[(S.1)] Remove a presolve record from the top of the stack.
\item[(S.2)] Transform the problem into the polyhedrally constrained setting.
\item[(S.3)] Undo the changes made to the model.
\item[(S.4)] Solve the optimization problem \eqnref{opt} using $\bar{x}$ to 
obtain $\bar{y}$.  The generated $(\bar{x}, \bar{y})$ solves the model before 
the presolve step was performed.
\item[(S.5)] Repeat until the presolve stack is empty.
\end{itemize}
}

\noindent
The optimization problem \eqnref{opt} is typically only in one dimension and 
is trivial to solve.  Care must be taken when calculating $N_X(\bar{x})$ 
because the algorithm may only find a solution to within a prespecified 
tolerance.  Therefore, variables within some tolerance of their bounds 
should be treated as if they are exactly on their bounds when constructing 
the normal cone.  The current default used in the implementation is the
convergence tolerance as specified by the algorithm.

The only case where two variables are involved in the optimization problem is 
when two inequalities are replaced by one equation.  The optimization problem 
in this case has an objective function equal to zero because at the solution
$A\bar{x} - b = 0$.  Therefore, a feasible point need only be generated.
In the presolved model, $0 \in F(x) - a\hat{y} + N_X(\bar{x})$
for the solution $(\bar{x}, \hat{y})$ given.  Without loss of generality, 
assume $Y_1 = \real_+$ and $Y_2 = \real_-$.  Select $y_1 \in Y_1 = \real_+$ 
and $y_2 \in Y_2 = \real_-$ such that $y_1 + y_2 = \hat{y}$.  These conditions
can always be trivially met.  Because the inclusion holds at $\hat{y}$, it
also holds for the $y_1$ and $y_2$ selected, which is then a feasible
point as required.

In the unfortunate case that the algorithm fails to solve the preprocessed
model, the optimization problem may have no solution either because it is 
infeasible or unbounded.  In this case, a value for the multiplier is chosen
in $Y$ such that the norm of the error in the box constrained representation 
is minimized given $\bar{x}$.  This greedy heuristic will lead to the best 
possible value in terms of the residual at each stage in the unrolling of the 
preprocessing steps, but not necessarily the least residual solution overall.

\section{Structural Implications}\label{sec:function}
The second phase of the preprocessor utilizes complementarity theory
to eliminate variables from the MCP.  The reductions documented are based 
on the rows and columns of the Jacobian, $F'$.  The main ingredient 
for the row-based rules is uniqueness.  If the value for a variable can be 
uniquely determined prior to solving the remainder of the problem, it is 
fixed and removed.  The column-based rules rely upon existence arguments.
Once a solution to the reduced model is known, a solution to the original
problem always exists and can be calculated.  Mechanisms developed include 
using interval evaluations, uncovering duplicate rows and columns, and 
exploiting special structure.  Note that when a row with zero elements and 
corresponding zero column are present in a model, the variable can
always be fixed at an appropriate value and removed.

\subsection{Intervals}
An interval evaluator determines the tightest possible $\underline{F}$ and 
$\bar{F}$ such that for all $x \in X$, $\underline{F} \le F(x) \le \bar{F}$.
For example, with a linear constraint, $F_i(x) = a^Tx - b$, the bounds
\begin{eqnarray*}
\underline{F}_i = \sum_{\{j \mid a_j > 0\}} a_j\ell_j + \sum_{\{j \mid a_j < 0\}} a_ju_j - b\\
\bar{F}_i = \sum_{\{j \mid a_j > 0\}} a_ju_j + \sum_{\{j \mid a_j < 0\}} a_j\ell_j - b\\
\end{eqnarray*}
can be used where $\ell_j$ and $u_j$ are the lower and upper bounds on variable 
$j$ respectively.  The range of a nonlinear function is dependent upon the
model and the bounds must be computed by a user supplied routine.

Using the ranges, variables in the model can be fixed.  If 
$\underline{F}_i > 0$, then $x_i$ must be fixed at its finite lower bound or 
the problem is infeasible.  Furthermore, if $\bar{F}_i < 0$, then $x_i$ must 
be fixed at its finite upper bound or the problem is infeasible.

Some of the constraints in the model will imply tighter bounds on the 
variables; i.e. a linear constraint.  These can be used by the interval 
evaluator to strengthen the ranges on other constraints, potentially
leading to more variables being fixed.  The tightened bounds are only 
used when calculating the intervals because modifying the bounds on the 
complementarity problem using this information could change the solution 
set of the problem.

\subsection{Duplicates}
Duplicate rows and columns can be very problematic for a solver.
By applying the same algorithm used in the polyhedrally constrained 
case (\secref{redundant}), two such linear rows or columns can be identified.
First consider the case of two duplicate rows in the problem.
Without loss of generality, the model can be written as:
\begin{eqnarray*}
0 \in \left[\begin{array}{c} F(x, y, z) \\ a^T(x, y, z) + b \\ a^T(x, y, z) + 
c\end{array}\right]
 + \left[\begin{array}{c}N_X(x) \\ N_Y(y) \\ N_Z(z) \end{array}\right].
\end{eqnarray*}
\tabref{duplicate} discusses the reductions that can be made.

\begin{table}
\centering
\caption{Duplicate Rows Cases}\label{tab:duplicate}
\begin{tabular}{lp{11cm}}
Case & Action \\
\hline
$Y = \real$ and $Z = \real$   & If $b \neq c$ the problem is infeasible.
                            Otherwise, nothing is done.\\
$Y = \real$ and $Z = \real_+$ & If $b > c$ fix $z$ at its lower bound.
                            Otherwise, if $b < c$, the problem is infeasible.\\
$Y = \real_+$ and $Z = \real_+$ & If $b > c$ fix $z$ at its lower bound.  
                            If $b < c$ fix $y$ at its lower bound.\\
$Y = \real_-$ and $Z = \real_+$ & If $b < c$, the problem is infeasible.  
                            Otherwise, nothing is done.\\
\end{tabular}
\end{table}

To remove column duplicates, one of the variables needs to be
free and the other must have two finite bounds.  The problem in this case
is:
\begin{eqnarray*}
0 \in F(x) + ay + az + N_{X \times \real \times [\ell, u]}(x, y, z),
\end{eqnarray*}
where $\ell$ and $u$ are the finite lower and upper bounds on $z$.
The reduction removes the $z$ variable from the problem and solves the 
reduced system to obtain $(\bar{x}, \bar{y})$.  
If $F_z(\bar{x}) + a_z\bar{y} > 0$, then $\hat{z} = \ell$.  Otherwise, 
$\hat{z} = u$.  Set $\hat{y} = \bar{y} - \hat{z}$.
Then $(\bar{x}, \hat{y}, \hat{z})$ solves the original problem 
as can be easily verified.

\subsection{Special Structure}
For a system of nonlinear equations, if the problem has the form:
\[
\left[\begin{array}{c} F(x) \\ G(x,y) \end{array}\right] = 
\left[\begin{array}{c} 0 \\ 0 \end{array}\right]
\]
with $F: \real^k \rightarrow \real^k$, then $F(x) = 0$ can be solved giving
$\bar{x}$ and then a $\bar{y}$ solving $G(\bar{x}, y) = 0$ can be found.
If $F(x) = 0$ has multiple solutions, this procedure may fail by finding an 
$\bar{x}$ for which $G(\bar{x},y) = 0$ has no solution.  For example,
consider $F(x) = x^2 - 1$ and $G(x,y) = x + y^2$.  $F(x) = 0$ has two solutions
$\bar{x} = 1$ and $\bar{x} = -1$.  Choosing $\bar{x} = 1$ leads to the
case where $G(\bar{x}, y)$ has no solution.
If $F(x) = 0$ has at most one solution, this case is precluded provided the 
original problem has a solution.  Similarly, the difficulty is alleviated if 
$G(\bar{x},y)$ has a solution $\bar{y}$ for all $\bar{x}$, since whatever 
$\bar{x}$ is found, the system $G(\bar{x}, y) = 0$ is solvable.  This 
section applies similar block reduction 
schemes to the mixed complementarity problem.

Consider a problem of the form:
\[
0 \in \left[\begin{array}{c} F(x) \\ G(x,y)\end{array}\right] +
\left[\begin{array}{c} N_X(x) \\ N_Y(y)\end{array}\right],
\]
where, as usual, $X$ and $Y$ are Cartesian products of intervals.
There are two sets of reductions that can be made.  If $0 \in F(x) + N_X(x)$ 
has a unique solution, $\bar{x}$, then $x$ can be fixed and the algorithm will
only work on the reduced problem $0 \in G(\bar{x}, y) + N_Y(y)$.  If $F(x)$ is
an affine function, i.e. $F(x) = Ax - b$, then it is known that 
$0 \in F(x) + N_X(x)$ has a unique solution if $A_X$, the normal map 
associated with this variational inequality, is coherently 
oriented \cite{robinson:normal}.  For example, when $X = \real_+^k$, this 
condition is equivalent to $A$ being a $P$-matrix.  For simple cases, coherent
orientation can be checked; e.g. when $k = 1$ or $2$.  In particular, when 
$F(x) = ax - b$ is a row singleton with a linear element on the diagonal, 
then coherent orientation is $a \neq 0$ when $X = \real$ and $a > 0$ in all 
other cases.  Satisfaction of this condition guarantees uniqueness of the
solution.  When $k = 2$, the condition is again that $A$ is a $P$-matrix, 
unless one or more of the intervals defining $X$ is $\real$.  In these later 
cases, the conditions are weaker.  \tabref{unique} summarizes all of the 
checks for coherent orientation.  The preprocessor identifies double blocks by 
finding a linear row with two elements, one of which is on the diagonal.  A 
check of the row corresponding to the other variable is performed to see if there 
is a doubleton block.

\begin{table}
\centering
\caption{Coherent Orientation Conditions}\label{tab:unique}
$X$ is assumed to be a Cartesian product of intervals with $\ell$ and $u$ 
being finite.

\begin{tabular}{ll}
Case & Coherent Orientation Condition \\
\hline
$X = \real$             & $A_{1,1} \ne 0$ \\
$X = [\ell,\infty)$      & $A_{1,1} > 0$ \\
$X = (-\infty,u]$     & $A_{1,1} > 0$ \\
$X = [\ell,u]$           & $A_{1,1} > 0$ \\
$X = \real \times \real$ & $det(A) \ne 0$ \\
$X = \real \times [\ell,\infty)$ & $det(A) \ne 0$ and $sign(det(A)) = sign(A_{1,1})$ \\
$X = \real \times (-\infty,u]$ & $det(A) \ne 0$ and $sign(det(A)) = sign(A_{1,1})$ \\
$X = \real \times [\ell, u]$ & $det(A) \ne 0$ and $sign(det(A)) = sign(A_{1,1})$ \\
All other cases & $det(A) > 0$ and $sign(A_{1,1}) = sign(A_{2,2}) = 1$
\end{tabular}
\end{table}

The other reduction to consider is where $0 \in G(x, y) + N_Y(y)$ has
a solution for all $x \in X$.  In this case,
$0 \in F(x) + N_X(x)$ is solved to find $\bar{x}$ and 
then a $\bar{y}$ satisfying $0 \in G(\bar{x}, y) + N_Y(y)$ is found.
Assume that $G(x, y)$ is linear in $y$, i.e.
$G(x,y) = H(x) + By$.  The coherent orientation conditions outlined above 
applied to $B$ suffice in this case as well, since they guarantee existence as
well as uniqueness.  However, to guarantee only existence, weaker conditions are
sufficient.
For $k = 1$ it is necessary and sufficient to have coherent
orientation or $Y$ compact.  When $k = 2$, the conditions are outlined in
\tabref{exist} and are derived from Theorem 2 in \cite{gowda:analysis} and
\cite{gowda:applications}.  The conditions given for the cases
where there is at least one free variable are necessary and sufficient
to guarantee existence for all $x \in X$.

\begin{table}
\centering
\caption{Existence Conditions}\label{tab:exist}
$Y$ is assumed to be a Cartesian product of intervals with $\ell$, 
$\tilde{\ell}$, $u$, and $\tilde{u}$ being finite.

\begin{tabular}{ll}
Case & Condition \\
\hline
$Y = \real$ & $B_{1,1} \neq 0$ \\
$Y = [\ell, \infty)$ & $B_{1,1} > 0$ \\
$Y = (-\infty, u]$ & $B_{1,1} > 0$ \\
$Y = [\ell, u]$ & nothing \\
$Y = \real \times \real$ & $det(B) \neq 0$ \\
$Y = \real \times [\ell, \infty)$ & $det(B) \neq 0$ and $sign(det(B)) = sign(B_{1,1})$ \\
$Y = \real \times (-\infty,u]$ & $det(B) \neq 0$ and $sign(det(B)) = sign(B_{1,1})$ \\
$Y = \real \times [\ell,u]$ & $B_{1,1} \neq 0$ \\
$Y = [\ell, \infty) \times [\tilde{\ell}, \infty)$ & $B > 0$ or \\
& $(det(B) > 0$ and $sign(B_{1,1}) = sign(B_{2,2}) = 1)$ \\
$Y = [\ell, \infty) \times (-\infty, \tilde{u}]$ & $B > 0$ or \\
& $(det(B) > 0$ and $sign(B_{1,1}) = sign(B_{2,2}) = 1)$ \\
$Y = [\ell, \infty) \times [\tilde{\ell},\tilde{u}]$ & $B_{1,1} \neq 0$ \\
$Y = (-\infty, u] \times (-\infty, \tilde{u}]$ & $B > 0$ or \\
& $(det(B) > 0$ and $sign(B_{1,1}) = sign(B_{2,2}) = 1)$ \\
$Y = (-\infty, u] \times [\tilde{\ell},\tilde{u}]$ & $B_{1,1} \neq 0$ \\
$Y = [\ell,u] \times [\tilde{\ell},\tilde{u}]$ & nothing
\end{tabular}
\end{table}

In the nonlinear setting, intervals on the Jacobian elements can
be used to verify conditions related to uniqueness and existence of
a solution.  For example in the single element case, if the value of
the Jacobian element is positive and uniformly bounded away from zero, i.e. it 
is a uniform $P$-function, then the existence and uniqueness is always 
guaranteed 
and the same substitutions can be performed.  Finding the solution becomes 
more difficult, as it involves solving a nonlinear problem.

\section{Computational Results}\label{sec:preproc-results}
The preprocessing algorithm implemented alternates between exploiting
the polyhedral structure and the functions.  Initially all possible 
reductions based on the polyhedral constraints are made.  Then all reductions 
based on the functional implications are made.  These two steps are
repeated until no changes are made to the model.

The potential for preprocessing is mainly limited to finding and
exploiting linear parts of the problem because interval evaluators are 
not available at present in the modeling language environments.  The
majority of the reductions made come from exploiting
polyhedral structure.  However, the reductions from
the second stage can also be significant to the success of the
algorithm.

The preprocessor was tested on three different sets of problems.  The
first test compares the performance of the MCP preprocessor to the 
one used by the commercial CPLEX code \cite{:cplex}.  Using the linear
programs contained in NETLIB \cite{gay:electronic}, the first order
conditions from linear programming were constructed and given to the 
MCP preprocessor.  CPLEX was given the original linear program.  
Reported in Tables \ref{tab:linear} and \ref{tab:linear1} are the sizes of 
the preprocessed models.  CPLEX is capable of performing aggregations, while
the MCP preprocessor currently does not.  Therefore, in the tables, the size
of the model produced by CPLEX both with aggregations (With) and without
aggregations (Without) are stated.  As evidenced by the table, the MCP 
preprocessor is competitive with CPLEX on linear programs when aggregations 
are not allowed.  One interesting point to note is that the problems
\texttt{fit*p} and \texttt{fit*d} are primal-dual pairs - the MCP preprocessor
generates an identical system in both cases.  Exploiting dual information 
in the \texttt{fit*p} problems significantly reduces the size of the
preprocessed models.

\begin{table}
\centering
\caption{Comparison of CPLEX and MCP preprocessor on NETLIB problems}\label{tab:linear}
\begin{tabular}{lcccc}
& & \multicolumn{2}{c}{CPLEX} \\
Model & Size & With & Without & MCP Preprocessor \\
\hline
adlittle & 153 & 147 & 147 & 146 \\
afiro & 59 & 48 & 52 & 56 \\
agg & 651 & 271 & 275 & 433 \\
agg2 & 818 & 530 & 538 & 743 \\
agg3 & 818 & 531 & 541 & 743 \\
bandm & 777 & 398 & 483 & 467 \\
beaconfd & 435 & 55 & 220 & 205 \\
blend & 157 & 108 & 140 & 149 \\
bnl1 & 1807 & 1443 & 1668 & 1670 \\
bnl2 & 5769 & 3031 & 4226 & 4341 \\
boeing1 & 909 & 711 & 713 & 720 \\
boeing2 & 320 & 281 & 281 & 292 \\
bore3d & 547 & 105 & 182 & 191 \\
brandy & 431 & 265 & 311 & 311 \\
capri & 608 & 383 & 547 & 547 \\
cycle & 4743 & 2700 & 3416 & 3884 \\
czprob & 4221 & 2904 & 3349 & 3430 \\
d2q06c & 7338 & 6450 & 6871 & 6286 \\
d6cube & 6588 & 5844 & 5867 & 6423 \\
degen2 & 978 & 855 & 977 & 974 \\
degen3 & 3321 & 3125 & 3321 & 3310 \\
dfl001 & 18301 & 13062 & 17091 & 16915 \\
e226 & 505 & 397 & 411 & 414 \\
etamacro & 1006 & 754 & 850 & 821 \\
fffff800 & 1378 & 933 & 983 & 1284 \\
finnis & 1066 & 739 & 786 & 808 \\
fit1d & 1050 & 1048 & 1048 & 1050 \\
fit1p & 2304 & 2054 & 2054 & 1050 \\
fit2d & 10525 & 10388 & 10388 & 10525 \\
fit2p & 16525 & 16525 & 16525 & 10525 \\
forplan & 554 & 466 & 476 & 483 \\
ganges & 2990 & 1202 & 2177 & 2466 \\
gfrd-pnc & 1708 & 1116 & 1656 & 1656 \\
greenbea & 7691 & 4055 & 5900 & 5763 \\
greenbeb & 7679 & 4044 & 5892 & 5738 \\
grow15 & 945 & 945 & 945 & 945 \\
grow22 & 1386 & 1386 & 1386 & 1386 \\
grow7 & 441 & 441 & 441 & 441 \\
israel & 316 & 304 & 304 & 304
\end{tabular}
\end{table}

\begin{table}
\centering
\caption{Comparison of CPLEX and MCP preprocessor on NETLIB problems (cont.)}\label{tab:linear1}
\begin{tabular}{lcccc}
& & \multicolumn{2}{c}{CPLEX} \\
Model & Size & With & Without & MCP Preprocessor \\
\hline
kb2 & 84 & 67 & 79 & 82 \\
lotfi & 461 & 399 & 399 & 408 \\
nesm & 3585 & 3325 & 3373 & 3440 \\
perold & 1937 & 1571 & 1769 & 1757 \\
pilot4 & 1380 & 1111 & 1200 & 1210 \\
sc105 & 207 & 117 & 207 & 207 \\
sc205 & 407 & 231 & 405 & 405 \\
sc50a & 97 & 57 & 97 & 97 \\
sc50b & 96 & 56 & 96 & 96 \\
scagr25 & 971 & 591 & 841 & 734 \\
scagr7 & 269 & 159 & 229 & 194 \\
scfxm1 & 787 & 612 & 694 & 698 \\
scfxm2 & 1574 & 1228 & 1388 & 1396 \\
scfxm3 & 2361 & 1844 & 2082 & 2094 \\
scorpion & 746 & 172 & 590 & 453 \\
scrs8 & 1659 & 913 & 1429 & 1438 \\
scsd1 & 837 & 837 & 837 & 817 \\
scsd6 & 1497 & 1497 & 1497 & 1481 \\
scsd8 & 3147 & 3147 & 3147 & 3135 \\
sctap1 & 780 & 608 & 608 & 660 \\
sctap2 & 2970 & 2303 & 2303 & 2500 \\
sctap3 & 3960 & 3111 & 3111 & 3340 \\
share1b & 342 & 297 & 315 & 310 \\
share2b & 175 & 168 & 172 & 172 \\
shell & 2061 & 1427 & 1935 & 1935 \\
ship04l & 2478 & 2174 & 2182 & 2208 \\
ship04s & 1818 & 1426 & 1482 & 1500 \\
ship08l & 4995 & 3569 & 3569 & 3611 \\
ship08s & 3099 & 1760 & 1858 & 1890 \\
ship12l & 6469 & 4756 & 4756 & 4790 \\
ship12s & 3805 & 2114 & 2258 & 2288 \\
stair & 741 & 512 & 740 & 740 \\
stocfor1 & 228 & 113 & 190 & 188 \\
stocfor2 & 4188 & 2474 & 3822 & 3825 \\
tuff & 878 & 514 & 738 & 788 \\
wood1p & 2838 & 1898 & 1898 & 1971 \\
\hline
Total & 181745 & 137202 & 155734 & 150753
\end{tabular}
\end{table}

A second test was performed using quadratic programming problems reformulated 
using the complementary slackness conditions \eqnref{qpcomp}.  Some artificial
quadratic programs were created for testing purposes from the NETLIB
collection.  A term of $\frac{1}{2}x^Tx$ was added to the objective function
and the resulting complementary slackness conditions were given to the
preprocessor and the PATH 4.x code.
\tabref{quadratic} reports the
size reductions and compares the solution times for PATH 4.x on the original 
and preprocessed models.  On the problems successfully preprocessed, the 
reductions in time are significant.  Some quadratic programs from
other sources were also tested.  In one of the models, \texttt{hwayoung}, 
over 70\% of the variables were removed by the preprocessor reducing the size 
from 46123 variables to 13655, leading to a significant reduction in
the total time needed to solve the problem.

\begin{table}
\centering
\caption{Comparison of PATH 4.x with and without preprocessing on QP models}\label{tab:quadratic}
\begin{tabular}{lcccc}
& \multicolumn{2}{c}{Original} & \multicolumn{2}{c}{Preprocessed}\\
Model & Size & Solution Time & Size & Solution Time  \\
\hline
agg & 651 & 6.6 & 454 & 1.0 \\
beaconfd & 435 & 1.1 & 283 & 0.7 \\
finnis & 1066 & 9.1 & 918 & 1.5 \\
lotfi & 461 & 5.9 & 434 & 0.9 \\
nesm & 3585 & 57.6 & 3481 & 53.0 \\
scorpion & 746 & 1.1 & 617 & 0.8 \\
ship08s & 3099 & 6.3 & 1966 & 3.3 \\
tuff & 878 & 4.3 & 849 & 3.9 \\
\end{tabular}
\end{table}

Finally, the models in MCPLIB \cite{dirkse.ferris:mcplib} were given
to the preprocessor.  The results on these models are less encouraging
than the other two tests because of a lack of linear problems
in the test set and the inability to obtain interval evaluations for
the nonlinear functions.  Many of the models did not benefit from
preprocessing.  However, some successes are reported in
\tabref{mixed}.  Note that the \texttt{explcp} model that is supposed
to display exponential behavior for Lemke's algorithm is completely
solved in the preprocessor.  The preprocessor for the \texttt{golanmcp} 
model removes 18 redundant rows.  The remaining problem solves without
any proximal perturbation, leading to the substantial 
reduction in solution time.

\begin{table}
\centering
\caption{Comparison of PATH 4.x with and without preprocessing on MCP models}\label{tab:mixed}
\begin{tabular}{lcccc}
& \multicolumn{2}{c}{Original} & \multicolumn{2}{c}{Preprocessed}\\
Model & Size & Solution Time & Size & Solution Time \\
\hline
electric & 158 & 1.3 & 149(143) & 0.5 (0.1) \\
explcp & 16 & 0.0 & 0 & 0.0 \\
forcebsm & 184 & 0.1 & 72 & 0.2 (0.1) \\
forcedsa & 186 & 0.1 & 70 & 0.1 (0.1) \\
golanmcp & 4321 & 80.9 & 4303 & 25.0 \\
merge & 9536 & 2254.6 & 8417 & 1954.2
\end{tabular}
\end{table}

In our initial testing, 
some of the preprocessing performed was
detrimental.  For example, the \texttt{force*} models became harder to solve
after preprocessing even though they were significantly reduced in
size.  We added an option that reports extra information to the
modeler regarding the orientation of equality constraints.  
This uncovered four models in MCPLIB that have equality constraints 
not formulated in a skew symmetric sense, namely the \texttt{electric}, 
\texttt{force*} and 
\texttt{lincont} models.  Rewriting to change this orientation gave 
the improved timings noted in parentheses in the table.

\section{Summary}
This chapter has developed a complete preprocessor for complementarity
problems and demonstrated the effectiveness of the procedure in reducing
the time spent by PATH 4.x to solve available test problems.  The preprocessor 
is currently available with both PATH 4.x and SEMI, which are the first codes 
for solving complementarity problems with preprocessing technology.  

We remark that an earlier version of this chapter appeared in 
\cite{ferris.munson:preprocessing}.

\chapter{Diagnostic Information}\label{chap:diagnostic}
Developing a practical model of a complex situation is a difficult task in 
which an approximate representation is initially constructed and then 
iteratively refined until an accurate formulation is obtained.  During the 
intermediate stages, the models generated have a tendency to be ill-defined, 
poorly conditioned, and/or singular.  Information generated by a solver can 
help the modeler to detect these problems, quickly locate the source, and make 
appropriate modifications to the model.  By preventing the propagation
of errors to successive models, the development cycle shortens and the final 
product becomes easier to solve with a more meaningful solution.

This chapter discusses statistics provided that can be used to find potential 
difficulties.
We begin by looking at the merit functions used to indicate if an iterate 
is close to the solution set of the given MCP in \secref{merit-functions}.  
We then utilize the merit function information in \secref{ill-defined} when 
we present the first of our problem areas, ill-defined models.  
\secref{poorly-scaled} addresses poorly-scaled models and the difficulties 
encountered in such circumstances.  Information provided to detect scaling 
problems are mentioned.  \secref{singular} concludes the trilogy of problems 
areas by looking at singular models.

We remark that the statistics generated are not a replacement for the 
modeler's knowledge about the application, rather it augments the available 
information to help them rapidly identify potential problem areas.

\section{Merit Functions}\label{sec:merit-functions}
A solver for complementarity problems typically employs a merit function to
indicate the closeness of the current iterate to the solution set.  The merit 
function is zero at a solution to the original problem and strictly positive 
otherwise.  Numerically, an algorithm terminates when the merit function is 
approximately equal to zero, thus possibly introducing spurious ``solutions''.

The modeler needs to be able to determine with some reasonable degree
of accuracy whether the algorithm terminated at solution or if it simply 
obtained a point satisfying the desired tolerances 
that is not close to the solution set.  For complementarity problems, we can 
provide several indicators with different characteristics to help make
such a determination.  If one of the indicators is not close to zero, then 
there is some evidence that the algorithm has not found a solution.  We note
that if all of the indicators are close to zero, we are reasonably sure we 
have found a solution.  However, the modeler has the final responsibility
to evaluate the ``solution'' and check that it makes sense for their 
application.

For the NCP, a standard merit function is 
\begin{eqnarray*}
\norm{(-x)_+, (-F(x))_+, [(x_i)_+(F_i(x))_+]_i}
\end{eqnarray*}
with the first two terms 
measuring the infeasibility of the current point and the last term indicating 
the complementarity error.  In this expression, we use $(\cdot)_+$ to 
represent the Euclidean projection of $x$ onto the nonnegative orthant, that 
is $(x)_+ = \max(x,0)$.  For the more general MCP, we can define a similar
function:
\[
  \norm{x-\pi(x),
        \left[\left(\frac{x_i - \ell_i}{\norm{\ell_i} + 1}\right)_+(F_i(x))_+\right]_i,
        \left[\left(\frac{u_i - x_i}{\norm{u_i} + 1}\right)_+(-F_i(x))_+\right]_i}
\]
where $\pi(x)$ represents the Euclidean projection of $x$ onto the set 
$[\ell,u]$.  We can see that if we have an NCP, the function is exactly 
the one previously given and for nonlinear systems of equations, this 
becomes $\norm{F(x)}$.

There are several reformulations of the MCP as systems of nonlinear (nonsmooth)
equations for which the corresponding residual is a natural merit function.
Some of these are as follows:
\begin{itemize}
\item Generalized Minimum Map: $x - \pi(x - F(x))$
\item Normal Map: $F(\pi(y)) + y - \pi(y)$
\item Fischer Function: $\Phi(x)$, where $\Phi_i(x) \define
\phi(x_i,F_i(x))$ with
\[
\phi(a,b) \define \sqrt{a+b} - a - b.
\]
Note that $\phi(a,b) = 0$ if and only if $0 \leq a \perp b \geq 0$.
\end{itemize}
In the context of nonlinear complementarity problems the generalized minimum
map corresponds to the classic minimum map $\min(x,F(x))$.  
Furthermore, for NCPs the 
minimum map and the Fischer function are both local error bounds and were
shown to be equivalent in \cite{tseng:growth}.  \figref{merit} in the 
subsequent section plots all of these merit functions for the ill-defined 
example discussed therein and highlights the differences between them.

The squared norm of $\Phi$, namely 
$\Psi(x) \define \frac{1}{2} \sum \phi(x_i,F_i)^2$, is continuously
differentiable on $\real^n$ provided $F$ itself is.  Therefore, the first order
optimality conditions for the unconstrained minimization of $\Psi(x)$, namely
$\nabla \Psi(x) = 0$ give another indication as to whether the point
under consideration is a solution of MCP. 

\section{Ill-Defined Models}\label{sec:ill-defined}
A problem can be ill-defined for several different reasons.  We concentrate
on the following particular cases.  We will call $F$ well-defined at 
$\bar{x} \in [\ell,u]$ if $\bar{x} \in \dom F$, the domain of $F$, and 
ill-defined at $\bar{x}$ otherwise.  Furthermore, we define $F$ to be 
well-defined near $\bar{x} \in [\ell,u]$ if there exists an open neighbourhood 
of $\bar{x}$, $\mathcal{N}(\bar{x})$, such that 
$[\ell,u]\cap\mathcal{N}(\bar{x}) \subseteq \dom F$.  By saying the 
function is well-defined near $\bar{x}$, we are simply stating that $F$ is 
defined for all $x \in [\ell,u]$ sufficiently close to $\bar{x}$.  A function
not well-defined near $\bar{x}$ is termed ill-defined near $\bar{x}$.

We will say that $F$ has a well-defined Jacobian at $\bar{x} \in [\ell,u]$ 
if there exists an open neighbourhood of $\bar{x}$, $\mathcal{N}(\bar{x})$, 
such that $\mathcal{N}(\bar{x}) \subseteq \dom F$ and $F$ is continuously 
differentiable on $\mathcal{N}(\bar{x})$.  Otherwise the function has an 
ill-defined Jacobian at $\bar{x}$.  We note that a well-defined Jacobian 
at $\bar{x}$ implies that the MCP has a well-defined function near 
$\bar{x}$, but the converse is not true.

Many solvers use both function and Jacobian information when attempting 
to solve an MCP.  Therefore, both of these definitions are relevant.  
We discuss cases where the function and Jacobian are ill-defined in the 
next two subsections.  We illustrate uses for the merit function information 
and final point statistics within the context of these problems.

\subsection{Function Undefined}
We begin with a one-dimensional problem for which $F$ is ill-defined at 
$x = 0$ as follows:
\[
\begin{array}{rcl}
  0 \leq x & \perp & \frac{1}{x} \geq 0.
\end{array}
\]
Here $x$ must be strictly positive because $\frac{1}{x}$ is undefined
at $x = 0$.  This condition implies that $F(x)$ must be equal to zero. Since 
$F(x)$ is strictly positive for all $x$ strictly positive, this problem has
no solution.

We are able to perform this analysis because the dimension of the problem 
is small.  In general, a modeler might not know a priori that a problem 
has no solution and might attempt to formulate and solve it.  For this 
particular example, we must specify an initial value for $x$, as the function
is undefined at the default starting point, $x = 0$.
The AMPL and GAMS interfaces would notice that $F$ is undefined at this
initial point, and terminate with an error message.
Note that such an error messages only indicates that the function $\frac{1}{x}$
is undefined as $x = 0$.  The complementarity problem might have a solution.

After setting a starting point and passing the model along to the
algorithm, we proceed to ``solve'' the model.  A portion 
of the output with relevant statistics about the solution is 
given in \figref{badf-output}.  At the end of the solve,
all of the merit 
functions given in \secref{merit-functions} are evaluated at the final point.  
The 
Normal Map merit function, and to a lesser extent, the complementarity 
error, indicate that the ``solution'' found does not necessarily solve the MCP.

\begin{figure}
\begin{verbatim}
FINAL STATISTICS
Inf-Norm of Complementarity . .  1.0000e+00 eqn: (F)
Inf-Norm of Normal Map. . . . .  1.1181e+16 eqn: (F)
Inf-Norm of Minimum Map . . . .  8.9441e-17 eqn: (F)
Inf-Norm of Fischer Function. .  8.9441e-17 eqn: (F)
Inf-Norm of Grad Fischer Fcn. .  8.9441e-17 eqn: (F)

FINAL POINT STATISTICS
Maximum of X. . . . . . . . . .  8.9441e-17 var: (X)
Maximum of F. . . . . . . . . .  1.1181e+16 eqn: (F)
Maximum of Grad F . . . . . . .  1.2501e+32 eqn: (F)
                                            var: (X)
\end{verbatim}
\caption{Output for Ill-Defined Function}
\label{fig:badf-output}
\end{figure}

To indicate the difference between the merit functions, \figref{merit} plots
them all for the simple example.  We note that as $x$ approaches positive 
infinity, numerically,  we are at a solution to the problem with respect to 
all of the merit functions except for the complementarity error, which 
remains equal to one.  As $x$ approaches zero, the merit functions diverge, 
also indicating that $x = 0$ is not a solution.  

\begin{figure}
\begin{center}
\includegraphics*[width=4in]{residual2.eps}
\end{center}
\caption{Merit Function Plot}
\label{fig:merit}
\end{figure}

The natural residual and Fischer function tend toward $0$ as $x \downarrow 0$.
From these measures, we might think $x = 0$ is the solution.  However,
as previously remarked $F$ is ill-defined at $x = 0$.  $F$ and $F'$ 
become very large, indicating that the function (and Jacobian) might not be 
well-defined.  We 
might be tempted to conclude that if one of the merit function indicators is 
not close to zero, then we have not found a solution.  This conclusion is 
not always the case.  When one of the indicators is non-zero, we have 
reservations about 
the solution, but we cannot eliminate the possibility that we are actually
close to a
solution.  If we slightly perturb the original problem to
\[
\begin{array}{rcl}
  0 \leq x & \perp & \frac{1}{x + \epsilon} \geq 0
\end{array}
\]
for a fixed $\epsilon > 0$, the function is well-defined over $\real^n_+$
and has a unique solution at $x = 0$.  In this case, by starting at 
$x > 0$ and sufficiently  small, all of the merit functions, with the 
exception of the Normal Map, indicate that we have solved the problem as is
shown by the output in \figref{goodf-output} for $\epsilon = 1*10^{-6}$ and
$x = 1*10^{-20}$.
\begin{figure}
\begin{verbatim}
FINAL STATISTICS
Inf-Norm of Complementarity . .  1.0000e-14 eqn: (G)
Inf-Norm of Normal Map. . . . .  1.0000e+06 eqn: (G)
Inf-Norm of Minimum Map . . . .  1.0000e-20 eqn: (G)
Inf-Norm of Fischer Function. .  1.0000e-20 eqn: (G)
Inf-Norm of Grad Fischer Fcn. .  1.0000e-20 eqn: (G)

FINAL POINT STATISTICS
Maximum of X. . . . . . . . . .  1.0000e-20 var: (X)
Maximum of F. . . . . . . . . .  1.0000e+06 eqn: (G)
Maximum of Grad F . . . . . . .  1.0000e+12 eqn: (G)
                                            var: (X)
\end{verbatim}
\caption{Output for Well-Defined Function}
\label{fig:goodf-output}
\end{figure}
In this case, the Normal Map is quite large and we might think 
that the function and Jacobian are undefined.
This example illustrates the point that all of these tests are not 
infallible.
The modeler still needs to do some detective work to determine if they
have found a solution or if the algorithm is converging to a point where the
function is ill-defined.

\subsection{Jacobian Undefined}
Since the algorithms use Newton-like method to solve the problems,
they also
need the Jacobian of $F$ to be well-defined.  One model for which the 
function is well-defined over $[\ell, u]$, but for which the Jacobian 
is undefined at the solution is: $0 \le x \perp -\sqrt{x} \ge 0$.  This 
model has a unique solution at $x = 0$.

Using the PATH 4.x algorithm and starting from the point $x = 1*10^{-14}$,
the output given in \figref{badJ-output} is generated.
\begin{figure}
\begin{verbatim}
FINAL STATISTICS
Inf-Norm of Complementarity . .  1.0000e-07 eqn: (F)
Inf-Norm of Normal Map. . . . .  1.0000e-07 eqn: (F)
Inf-Norm of Minimum Map . . . .  1.0000e-07 eqn: (F)
Inf-Norm of Fischer Function. .  2.0000e-07 eqn: (F)
Inf-Norm of Grad FB Function. .  2.0000e+00 eqn: (F)

FINAL POINT STATISTICS
Maximum of X. . . . . . . . . .  1.0000e-14 var: (X)
Maximum of F. . . . . . . . . .  1.0000e-07 eqn: (F)
Maximum of Grad F . . . . . . .  5.0000e+06 eqn: (F)
                                            var: (X)
\end{verbatim}
\caption{Output for Ill-Defined Jacobian}
\label{fig:badJ-output}
\end{figure}
We can see the that gradient of the Fischer Function is nonzero and the 
Jacobian is beginning to become large.  These conditions indicate that the 
Jacobian is undefined at the solution.
It is therefore important for a modeler to inspect the given
output to guard against such problems.

If we start from $x = 0$, the algorithms correctly report that we are 
at the solution.  Even though the entries in the Jacobian are undefined 
at this point, some of the interfaces will not return a domain violation.
This problem with the Jacobian is therefore undetectable by the algorithm.

\section{Poorly Scaled Models}
\label{sec:poorly-scaled}
Problems which are well-defined can have various numerical problems that
can impede the algorithm's convergence.  One particular problem is a badly
scaled Jacobian.  In such cases, we can obtain a poor ``Newton'' direction
because of numerical problems introduced in the linear algebra performed.
This problem can also lead the code to a point from which it cannot recover.  

The final model given to the solver should be scaled such that we avoid 
numerical difficulties in the linear algebra.  The statistics provided
can be used to iteratively refine the model so that we eventually end up with 
a well-scaled problem.  We note that we only calculate our scaling statistics
at the starting point provided.  For nonlinear problems these statistics
may not be indicative of the overall scaling of the model.  Model specific 
knowledge is very important when we have a nonlinear problem because it can 
be used to appropriately scale the model to achieve a desired result.

We look at the \texttt{titan} model in MCPLIB, that has 
some scaling problems.  The relevant output for the original code
is given in \figref{poor-scaled}.
\begin{figure}
\begin{verbatim}
INITIAL POINT STATISTICS
Maximum of X. . . . . . . . . .  4.1279e+06 var: (w.29)
Maximum of F. . . . . . . . . .  2.2516e+00 eqn: (a1.33)
Maximum of Grad F . . . . . . .  6.7753e+06 eqn: (a1.29)
                                            var: (x1.29)

INITIAL JACOBIAN NORM STATISTICS
Maximum Row Norm. . . . . . . .  9.4504e+06 eqn: (a2.29)
Minimum Row Norm. . . . . . . .  2.7680e-03 eqn: (g.10)
Maximum Column Norm . . . . . .  9.4504e+06 var: (x2.29)
Minimum Column Norm . . . . . .  1.3840e-03 var: (w.10)
\end{verbatim}
\caption{Output - Poorly Scaled Model}
\label{fig:poor-scaled}
\end{figure}
The maximum row norm is defined as 
\[ \max_{1 \le i \le n} \sum_{1 \le j \le n} \mid [F'(x)]_{ij}
\mid \] 
and the minimum row norm is 
\[ \min_{1 \le i \le n} \sum_{1 \le j \le n} \mid [F'(x)]_{ij}
\mid . \]
Similar definitions are used for the column norm.
The norm numbers for this particular example are not extremely large,
but we can nevertheless improve the scaling.  We first decided to reduce the
magnitude of the \texttt{a2} block of equations as indicated by the output.   
After scaling the \texttt{a2} block, we re-ran the code
found additional blocks of equations and variables that also needed scaling. 
After scaling all of these blocks of equations in the model, we have improved 
the scaling statistics which are given in \figref{well-scaled} for the
new model.
\begin{figure}
\begin{verbatim}
INITIAL POINT STATISTICS
Maximum of X. . . . . . . . . .  1.0750e+03 var: (x1.49)
Maximum of F. . . . . . . . . .  3.9829e-01 eqn: (g.10)
Maximum of Grad F . . . . . . .  6.7753e+03 eqn: (a1.29)
                                            var: (x1.29)

INITIAL JACOBIAN NORM STATISTICS
Maximum Row Norm. . . . . . . .  9.4524e+03 eqn: (a2.29)
Minimum Row Norm. . . . . . . .  2.7680e+00 eqn: (g.10)
Maximum Column Norm . . . . . .  9.4904e+03 var: (x2.29)
Minimum Column Norm . . . . . .  1.3840e-01 var: (w.10)
\end{verbatim}
\caption{Output - Well-Scaled Model}
\label{fig:well-scaled}
\end{figure}
For this particular problem the PATH algorithm could not solve the unscaled 
model, while it can find a solution to the scaled model.

Not all of the numerical problems are directly attributable to poorly
scaled models.  Problems for which the Jacobian of the active constraints
is singular or nearly singular can also cause numerical difficulty as
illustrated next.

\section{Singular Models}
\label{sec:singular}

Assuming that the problem is well-defined and properly scaled, we can
still have a Jacobian for which the active constraints are singular or
nearly singular (i.e. it is ill-conditioned).  When problems are singular or 
nearly singular, we are also likely to have numerical problems.  As a result 
the ``Newton'' direction obtained from the linear problem solver can be 
very bad.  Various heuristics are used by algorithms that
attempt to remove the singularity problems from the model.
However, it is most often beneficial for solver robustness to remove
singularities if possible.

The easiest problems to detect are those for which the Jacobian has
zero rows and columns.  A simple problem for which we have zero rows
and columns is:
\[
  \begin{array}{rcl}
  -2 \leq x \leq 2 & \perp & -x^2 + 1.
  \end{array}
\]
Note that the Jacobian, $-2x$, is non-singular at all three solutions, but
singular at the point $x = 0$.  Output for this model starting at 
$x = 0$ is given in \figref{zero}.
\begin{figure}
\begin{verbatim}
INITIAL POINT STATISTICS
Zero column of order. . . . . .  0.0000e+00 var: (X)
Zero row of order . . . . . . .  0.0000e+00 eqn: (F)
Total zero columns. . . . . . .  1
Total zero rows . . . . . . . .  1
Maximum of F. . . . . . . . . .  1.0000e+00 eqn: (F)
Maximum of Grad F . . . . . . .  0.0000e+00 eqn: (F)
                                            var: (X)
\end{verbatim}
\caption{Output - Zero Rows and Columns}
\label{fig:zero}
\end{figure}
We display in the code the variables and equations for which the
row/column in the Jacobian is close to zero.  These situations are 
problematic and for nonlinear problems likely stem from the modeler providing 
an inappropriate starting point or fixing some variables resulting in some
equations becoming constant.  We note that the solver may
perform well in the presence of zero rows and/or columns, but the 
modeler should make sure that these are what was intended.

Singularities in the model can also be detected by the linear solver.
This in itself is a hard problem and prone to error.  For matrices
which are poorly scaled, we can incorrectly identify ``linearly dependent''
rows because of numerical problems.  Typically, singularity does not cause 
a lot of problems and the algorithm can handle the situation appropriately.
However, an excessive number of singularities are cause for concern.
A further indication of possible singularities at the solution is the lack of
quadratic convergence to the solution.

\section{Summary}
This chapter provided documentation on diagnostic information used to identify
ill-defined, poorly conditioned, and/or singular models.
We remark that an earlier version of this chapter appeared in 
\cite{ferris.munson:case}.

