
% MODIFY THE FOLLOWING TO INCLUDE THE LECTURE #, TOPIC, DATE, and YOUR NAME
\lecturetitle{9}{Linear Programming 2}{26 Sep 97}{Steven Czerwinski}

% INCLUDE AN OUTLINE OF THE CLASS TOPICS
\begin{itemize} 
  \item Duality

  \item Ellipsoid Algorithm
  \item Interior Point Methods
    \begin{enumerate}
    \item Affine scaling
    \item Potential reduction
    \item Central trajectory
    \end{enumerate}
\end{itemize}

\section{Duality}
\label{sec:duality}

{\em Duality} refers to the property that linear programming problems
can be stated both as minimization and maximization problems.  For any
problem, there is a {\em primal} and a {\em dual}.  
The primal usually refers to the most natural way we describe the original
problem.  This can either be a minimization or maximization problem,
depending on the specifics.  
The dual represents an alternative way to specify the original
problem, such that it is a minimization problem if the primal is a
maximization and vice versa.  Solving the dual is equivalent to
solving the original problem.

For example, suppose we have the following:
\begin{quote}
\begin{tabular}{ll}
{\bf Primal problem:} & Maximize $z=cx$ \\
& such that $Ax \leq b$ and $x \geq 0$ \\
& ($n$ equations, $m$ variables) \\
\end{tabular}
\end{quote}
Then, for our dual, we have: ({\em notice the change in $m$ and $n$})
\begin{quote}
\begin{tabular}{ll}
{\bf Dual problem:} & Minimize $z=yb$ \\
& such that $yA \geq c$ and $y \geq 0$ \\
& ($m$ equations, $n$ variables) \\
\end{tabular}
\end{quote}

So why do we care about the dual for a problem?  One reason is that we often
find that the dual is easier to solve.  
The dual can also be used for {\em sensitivity analysis} which
determines how much the solution would change with changes in $c,A$ or $b$.
Finally, it can be used in conjunction with the primal problem to help solve
some of the interior point methods.

%The Dual simplex method also
%uses duality to pick which problem will be quicker for simplex to
%solve.  Since simplex is dependent on the structure of the matrix it
%is using to solve the problem, a sensitivity analysis can show which
%version of the problem is better for simplex.

% Okay, but what about this ellipsoid alg bit

\subsection{Duality: General Case}

In converting maximization problems to minimization problems and vice
versa, we find that every variable in the primal produces a constraint
in the dual.  Likewise, every constraint in the primal produces a
variable in the dual.  

We also find there is correspondence between the variable's
restrictions in the primal, with the constraint's equality sign in the
dual.  For example, if a variable must be greater than 0 in a
minimization problem, then the constraint the variable produces in the dual 
will be stated as a ``less than'' relationship.  Similarly, there is
a correspondence between the constraints' equality signs in the primal
to the restrictions on the variables in the dual.
Table~\ref{tbl:dual-convert} summaries these correspondences.  

\begin{table}
\centering
\begin{tabular}{|ccc|} \hline
{\bf Maximization Problem} & & {\bf Minimization Problem} \\ \hline
Constraints & & Variables \\
$\leq$ & $\Leftrightarrow$ & $\geq 0$ \\
$\geq$ & $\Leftrightarrow$ & $\leq 0$ \\
$=$ & $\Leftrightarrow$ & unrestricted \\
& &\\
Variables & & Constraints \\
$\geq 0$ & $\Leftrightarrow$ & $\geq$ \\
$\leq 0$ & $\Leftrightarrow$ & $\leq$ \\
unrestricted & $\Leftrightarrow$ & $=$ \\ \hline
\end{tabular}
\caption{Rules for converting between Primals and Duals}
\label{tbl:dual-convert}
\end{table}

\subsection{Example}

We are given the following maximization problem:
\begin{quote}
\begin{tabular}{rl}
maximize: & $4x_{1} + 2x_{2} - x_{3}$ \\
& \\
s.t. & $ x_{1} + x_{2} + x_{3} = 20 $\\
     & $ 2x_{1} - x_{2} \geq 6 $\\
     & $ 3x_{1} + 2x_{2} + x_{3} \leq 40 $\\
     & $ x_{1}, x_{2} \geq 0 $\\
\end{tabular}
\end{quote}
From this description, we easily determine the following:

\begin{eqnarray*}
x &=& \left[ \begin{array}{l} x_{1} \\ x_{2} \\ x_{3} \end{array}
\right] \\
c &=& \left[ \begin{array}{lll} 4 & 2& -1 \end{array} \right] \\
b &=&  \left[ \begin{array}{l} 20 \\ 6 \\ 40 \end{array} \right] \\
A &=& \left[ \begin{array}{lll}
1 & 1 & 1 \\
2 & -1 & 0 \\
3 & 2 & 1 \end{array}
\right]
\end{eqnarray*}  

We can now produce the dual:
\begin{quote}
\begin{tabular}{rl}
minimize: & $20Y_{1} + 6Y_{2} - 40Y_{3}$ \\
& \\
s.t. & $ Y_{1} + 2Y_{2} + 3Y_{3} \geq 4 $\\
     & $ Y_{1} - Y_{2} + 2Y_{3} \geq 2 $\\
     & $ Y_{1}+ Y_{3} = -1 $\\
     & $ Y_{2} \leq 0, Y_{3} \geq 0 $\\
\end{tabular}
\end{quote}
The coefficients come from transposing $A$.  To determine the
equality/inequality type we use the rules in
Table~\ref{tbl:dual-convert}.  For example the first inequality
$ Y_{1} + 2Y_{2} + 3Y_{3} \geq 4 $ corresponds to variable $x_1$ in
the primal.  Since $x_1 \geq 0$ in the primal, the inequality type
is also $\geq$ in the dual.

Notice, since the first constraint is produce by multiplying $A$'s
first column vector with $Y$, it corresponds to the $x_{1}$ variable
in the primal (because the first column in $A$ contains all of
$x_{1}$'s coefficients).  Therefore, the first constraint in the dual
must have an ``equal'' relationship since $x_{1}$ is unrestricted in
the primal.

\subsection{Duality Theorem}

The Duality theorem states:

If $\overline{x}$ is feasible for {\bf P} (given below) and
$\overline{y}$ is feasible for {\bf D} then: $c\overline{x} \leq
\overline{y}b$ and at optimality, $c\overline{x} = \overline{y}b$.  We
call $\overline{y}b - c\overline{x}$ the {\em Duality Gap}.
The magnitude of the duality gap is a good way to determine how close
we are to an optimal solution and is used in some of the interior
point methods.

\begin{quote}
\begin{tabular}{rl}
{\bf P}:  Max& $z = cx $ \\
& \\
s.t. & $Ax \leq b$ \\
     & $x \geq 0$ \\
& \\
{\bf D}:  Min & $Z = yb$ \\
s.t. & $yA \geq c$ \\
     & $y \geq 0$
\end{tabular}
\end{quote}

\section{Ellipsoid Algorithm}
\index{Ellipsoid Algorithm: Khachian}

Created by Khachian in '79, the ellipsoid algorithm was the first
known polynomial-time algorithm for linear programming.  (There were
other poly-time algorithms around, but they had not been proven to be
so.)
The ellipsoid algorithm does not solve linear programming problems of
the general form, but rather it solves any problem stated in
the following way:  

Find $x$ subject to $Ax \leq b$ (i.e. find any feasible
solution).  

However, we will show that
any linear programming problem can be reduced to this form.

If $L$ equals the number of bits used to represent $A$ and $b$, then
the running time of the algorithm is bound by $O(n^{4}L)$.
Unfortunately the algorithm takes this much time for most problems (it
is very unusual that it finds a solution any faster).  Simplex's worse
case running time is larger than $O(n^{4}L)$, but for most problems it
does much better.  Therefore Simplex ends up being better than the Ellipsoid
algorithm in practice for most problems even though its worst-case
running time is worse.

\subsection{Description of algorithm}

The idea behind the ellipsoid algorithm is that we start with some
ellipsoid $E_{0}$ which encloses an area in the solution space, and is
guaranteed to contain the feasible region of the solution.  Basically,
at each iteration of the algorithm, we create new ellipsoids $E_{k}$
that are smaller than the previous, but still are guaranteed to 
contain the feasible region.  Eventually, we will find a point that
satisfies the problem (recall that we are satisfied with any feasible
point).

Specifically, for each iteration $k$ of the algorithm we do the
following:
\begin{enumerate}
\item Determine $c_{k} = $ the center of $E_{k}$.
\item Determine $a_{k} = $ the constraint most violated by $c_{k}$.
\item Use the $a_{k}$ hyperplane to split $E_{k}$ into two regions.
One of these regions $R_1$ contains $c_{k}$, and the other $R_2$
contains the feasible region.
\item We then construct $E_{k+1}$ such that it contains 
all points in $E_k \cap R_2$, but is smaller than $E_k$.
\end{enumerate}
Figure~\ref{fig:ellp-alg} illustrates these steps.

\begin{figure}
\centerline{\PSbox{../lecture09/ellip-alg2.eps}{5.1in}{3.75in}}
\caption{An iteration in the Ellipsoid Algorithm}
\label{fig:ellp-alg}
\end{figure}

For $n$ variables
Khachian showed the following property of this algorithm:

\[ \frac{Volume(E_{k+1})}{Volume(E_{k})} =
\frac{1}{2^{\frac{1}{2n+1}}} \]

\subsection{Reduction from general case}

As stated, we can convert the general form of linear programming
problems to a form suitable to use the ellipsoid algorithm on.

Suppose we have the following problem stated in the general form:
Maximize $cx$ subject to $Ax \leq b$, and $x \geq 0$.  We can convert
it to the form $A^{\prime} x^{\prime} \leq b^{\prime} $ as shown
below:

Find any $x$, $y$ such that
\begin{eqnarray}
Ax & \leq & b  \label{eq:e1}\\
-x & \leq & 0  \label{eq:e2}\\
yA & \leq & -c \label{eq:e3}\\
-y & \leq & 0 \label{eq:e4}\\
-cx + by & \leq & 0 \label{eq:e5}
\end{eqnarray}
Here $x$ and $y$ together form $x^{\prime}$, the combination of the
left sides of the equations form $A^{\prime}$, and the combination of
the right sides of equations form $b^{\prime}$.  Equations
(\ref{eq:e1}) and (\ref{eq:e2}) come from the original problem, and
equations (\ref{eq:e3}) and (\ref{eq:e4}) come from its dual.
Equation (\ref{eq:e5}) will only be true at a single point when the original
problem is optimized (the duality gap must be non-negative).  In this
way, if any solution to $A^{\prime} x^{\prime}
\leq b^{\prime} $ is found, then it will be the optimal solution for
the original problem.

\section{Interior Point Methods}

Interior point methods are considered the first practical polynomial
time methods for linear programming.  Unlike
simplex, the interior point methods do not touch any boundaries until
the end when they find the optimal solution.
Instead,
these methods work by traveling through the interior of the feasible
region (see Figure~\ref{fig:int-point}) by using a combination of the
following two terms:

\begin{description}
\item[Optimizing term]  This moves the algorithm towards the
objective.
\item[Centering term]  This keeps it away from the boundaries.
\end{description}

\noindent
Generally, interior point methods take many fewer steps than the
simplex method (often less than 100), but each step is considerably
more expensive.  The overall relative cost of simplex and interior
point methods is very problem specific and is effected by properties
beyond just the size of the problem.
For example the times depend strongly on the matrix structure.

\begin{figure}
\centerline{\PSbox{../lecture09/int-point.eps}{3.2in}{2in}}
\caption{Walking on the interior points}
\label{fig:int-point}
\end{figure}

Interior point methods have been used since the 50s for nonlinear
programming.  Roughly categorized, the methods include the following:
\begin{quote}
\begin{tabular}{l|l}
{\bf Name} & {\bf Number of iterations} \\ \hline
Affine scaling & no known time bounds \\
Potential reduction & $O(nL)$ \\
Central trajectory & $O(\sqrt{x}L)$ \\
\end{tabular}
\end{quote}
Karmarkar proved polynomial time bounds for his method (potential
reduction) in 1984.  Table~\ref{tbl:int-times} shows some example
running times for a central trajectory method (Lustig, Marsten,
and Shanno 94) on four different benchmark problems.  As the table shows,
the time spent in each iteration is greatly influenced by the number
of Cholesky non-zeros present in the matrix.

\begin{table}
\centering
\begin{tabular}{|r|cccc|}
\hline
& {\bf Fuel} & {\bf Continent} & {\bf Car} & {\bf Initial} \\ \hline
{\bf Size} & 13x31k & 9x57k& 43x107k & 19x12k \\
{\bf Non-zeros} & 186k & 189k & 183k & 80k\\
{\bf Iterations} & 66& 64& 53 & 58\\
{\bf Cholesky non-zeros} & 1.2M & .3M & .2M & 6.7M\\
{\bf Time} & 2364s & 771s & 645s & 9252s\\
\hline
\end{tabular}
\caption{Performance of the central trajectory method on 4 benchmarks}
\label{tbl:int-times}
\end{table}

In the following discussion we will assume we are trying to solve a
linear program of the form 
\begin{quote}
\begin{tabular}{rl}
minimize: & $c x$ \\
subject to: & $Ax = b$\\
            & $x \geq 0$
\end{tabular}
\end{quote}

\subsection{Centering}

There are several different techniques to compute the centering term
for interior point methods.  This section describes a few of those
techniques.

\subsubsection{Analytical center}

We can use analytical centering by minimizing 
\[ \sum_{i=1}^{m} \lg x_{i} \]
If $x_{i} = 0$ are the boundaries of the feasible region (we assume
they are), then $x_{i}$ are the distances to the boundaries from our
current point.  See Figure~\ref{fig:ana-bound}. 


\begin{figure}
\centerline{\PSbox{../lecture09/ana-bound.eps}{3.5in}{3in}}
\caption{The analytical center}
\label{fig:ana-bound}
\end{figure}

\subsubsection{Elliptical Scaling}

In elliptical scaling, we start with a small ellipsoid centered around
the current point.  We then expand the ellipsoid in each dimension of the
problem.  The expansion in a dimension stops when the ellipsoid
touches a constraint.  Since our inequality constraints
are of the form $x \geq 0$, the equation will be of the form:

\[ \frac{x_{1}^{2}}{c_{1}^{2}} + \frac{x_{2}^{2}}{c_{2}^{2}} = 1 \]

The ellipsoid will then be a measure of how uncentered the point is
and can be used to scale space to bias our selection of a direction
toward the longer dimensions of the ellipsoid, and thus
away from constraints which are close by (the shorter dimensions).
See Figure~\ref{fig:elp-scale}.
We will see how this is used in the affine scaling method.


\begin{figure}
\centerline{\PSbox{../lecture09/elp-scale.eps}{2in}{2.3in}}
\caption{Elliptical Scaling}
\label{fig:elp-scale}
\end{figure}

\subsection{Optimizing techniques}

\subsubsection{Projection}
\label{section:projection}

If our problem is constrained to $Ax = b$, then we know that our
feasible region is contained in a hyperplane.  If $A$ has $n$ rows and
$m$ columns, then it is contained in a $m - n$ dimensional hyperplane
in $m$ dimensional space. (Assuming the rows in $A$ are independent.)
The hyperplane basis is called the {\em null space} of $A$
We can also view $A$ as being the ``slopes'' for the hyperplane, while
$b$ represents the ``offsets''.

We can use this to find the optimal solution for the problem.   If
$\vec{c}$ is the direction of the optimal solution, then we should
take a step in the direction of $\vec{d}$.  See Figure~\ref{fig:proj}.
We must do this because $\vec{c}$ would take us off the hyperplane and
therefore, out of the feasible region.  Instead, we project $c$ into
the hyperplane, and step in that direction.

Here's how we compute $\vec{d}$:

\begin{eqnarray*}
\vec{d} & = &\vec{c} - \vec{n} \\
& = & c - A^{{\tt T}}(A A^{\tt{T}})^{-1}Ac \\
& = & (I -  A^{\tt{T}}(A A^{\tt{T}})^{-1}A)c \\
& = & P \cdot c
\end{eqnarray*}    
where $P = (I -  A^{\tt{T}}(A A^{\tt{T}})^{-1}A)$, the projection
matrix.

\begin{figure}
\centerline{\PSbox{../lecture09/proj.eps}{3in}{1.5in}}
\caption{Projection into the null space}
\label{fig:proj}
\end{figure}

Calculating $P \cdot c$ (the projection) must be done at each
iteration of the method since, as we will see, $A$ will change.  So,
solving for the projection matrix is a major computational concern.

We find that it is more efficient to solve $AA^{\tt{T}}w = Ac$ for
$w$.   We can then use $w$ to calculate the projection because of the
following relationship:

\[ Pc = (I -  A^{\tt{T}}(A A^{\tt{T}})^{-1}A)c = c -
A^{\tt{T}}w \]

Computing $AA^{\tt{T}}w = Ac$ is more efficient because it can be
solved using a sparse solver (e.g. Cholesky Factorization).
Otherwise, calculating $P$ directly can generate a dense matrix,
which is problematic both in terms of runtime and in the
space it would require.

\subsubsection{Newton's Method}

We can use Newton's Method to help determine the direction in which we
should travel in the interior of the feasible region to optimize the
solution.  Often, the function we are trying to optimize will not be
linear in $x$, and could be composed of transcendental functions
such as $\lg or \ln$, as in the central trajectory case.  We can
approximate the minimum of these functions using a Newton's method, and more easily
determine the direction of the optimal solution.  We can see how this
is done through a little math.

First, we will just look at the general case of minimizing some
1-dimensional function $f(x)$.  We can write $f(x)$ as a Taylor
expansion around some point $x_{0}$.

\[ f(x) \approx f(x_{0}) + f^{\prime} (x_{0})(x - x_{0}) + \frac{1}{2} 
f^{\prime \prime} (x_{0})(x - x_{0})^{2} \]

To get the minimum of the function, we take the derivative and set it
to 0, which gives:

\begin{eqnarray*}
0 & = & f^{\prime} (x_{0}) + f^{\prime \prime}(x_{0})(x - x_{0}) \\
x & = & x_{0} - \frac{ f^{\prime} (x_{0})}{ f^{\prime \prime}(x_{0})}
\end{eqnarray*}

In multiple dimensions we can write this in matrix form as

\[ x - x_c = \vec{d} = (F(x_c))^{-1} \nabla
f(x_c)^{\tt{T}} \]
where $f(x)$ is now a function of $n$ variables, $x_c$ is
the point we are expanding about, and $F(x_c) =
\nabla \times \nabla f(x)$ evaluated at $x_c$,
an $n\times n$ matrix which is called the Hessian matrix, i.e.
\[F(x) = \left[
\begin{array}{cccc}
\frac{\partial f}{\partial x_1 x_1} &
\frac{\partial f}{\partial x_1 x_2} &
\cdots &
\frac{\partial f}{\partial x_1 x_n} \\
\frac{\partial f}{\partial x_2 x_1} &
\frac{\partial f}{\partial x_2 x_2} &
\cdots &
\frac{\partial f}{\partial x_2 x_n} \\
& & \cdots & \\
\frac{\partial f}{\partial x_n x_1} &
\frac{\partial f}{\partial x_n x_2} &
\cdots &
\frac{\partial f}{\partial x_n x_n} 
\end{array}
\right] \]
$\nabla f(x_c)$ is a vector of the first
derivatives of $f(x)$ evaluated at $x_c$.

If we constrain the newton step to be on a hyperplane (i.e. $Ax = b$),
it is called a {\em constrained Newton step}.
In general we first pick a direction $\vec{d}$ based on a Newton
step (or other method) and then project the step onto the
constraining hyperplane, as discussed in the previous section.


\subsection{Affine scaling method}
\index{Interior Point: Affine Scaling}

The affine scaling method is basically a biased steepest descent
search (biased away from the boundaries of the feasible region).
Suppose we had the following problem to solve:
\begin{quote}
\begin{tabular}{rl}
minimize: & $cx$ \\
s.t. & $Ax = b$, \\
& $x \geq 0$
\end{tabular}
\end{quote}

To solve this using the affine scaling method, at each
iteration we solve:
\begin{quote}
\begin{tabular}{rl}
minimize: & $cz$ \\
s.t. & $Az = 0$, \\
& $ zD^{-2}z \leq 1$
\end{tabular}
\end{quote}
where $D$ is a matrix whose only non-zero terms are along its main
diagonal.  The terms on this diagonal are the coordinates of the
current point: $x_{1}, x_{2}, x_{3}, \cdots$.  The equation $ zD^{-2}z
\leq 1$ is the definition for the interior of an ellipsoid, as
described earlier (see Figure~\ref{fig:elp-scale}).

We update our $x$ at each iteration by setting $x_{k+1} = x_{k} +
\alpha * z $, where $z$ is the direction we found in the above
minimization problem, and $\alpha$ is selected to bring us close
to the boundary, but not all the way there.

In this process, the $cz$ minimization acts as the optimization term,
while the $ zD^{-2}z \leq 1$ acts as the centering term.  In order to
understand how this works, we should realize $Ax = b$ is a slice of
the ellipsoid.  In fact, it is the same hyperplane that we described
in the ellipsoid scaling method.  In order to minimize the function,
we should step in the $c^{\prime}$ direction, where $c^{\prime}$ is
the projection of $c$ into the hyperplane.  However, the affine
scaling method steps in the direction of $z$.  

How do we get $z$?  We construct hyperplanes perpendicular to $c$, and
keep moving them back until we find one such hyperplane
that intersects the ellipsoid only once.  We then construct $z$ to be
the difference between the current point and the intersection point.
By stepping in the direction of $z$, we are centering the point
because $z$ will tend to point to constraints that are far away.  In
fact, as the current point comes closer and closer to a constrain,
$z$ will point further and further away from that constraint.

See Figure~\ref{fig:affine} for a picture of this construction of $z$,
$c^{\prime}$, and the hyperplanes.

\begin{figure}
\centerline{\PSbox{../lecture09/affine.eps}{3.1in}{3.5in}}
\caption{An iteration of the Affine Scaling Method}
\label{fig:affine}
\end{figure}

\subsubsection{Computation at each iteration}

In order to ease the computation of $z$, we make a substitution of
variables, by setting $z = Dz^{\prime}$.  Then, we have
\begin{quote}
\begin{tabular}{rl}
minimize: & $cDz^{\prime}$ \\
s.t. & $ADz^{\prime} = 0$, \\
& $ z^{\prime \tt{T}}D^{\tt{T}}D^{-2}Dz^{\prime} \leq 1$
\end{tabular}
\end{quote}
The last equation simplifies to $zz \leq 1$, which is just a
sphere.  With this in mind, we can say our direction $z^{\prime}$ is
unbiased and throw out the $zz \leq 1$ constraint.
To find the minimum we now just take a steepest descent step on $cDz^{\prime}$
that keeps us within the hyperplane defined by $ADz^{\prime} = 0$.
This leads us to,

\[ z^{\prime} = - \underbrace{(I -
B^{\tt{T}}(BB^{\tt{T}})^{-1}B)\cdot cD}_{\tt{projection}} \]
where $B = AD$.  To get $z$ in our original space we just
substitute back $z = Dz'$.   The vector $z'$ is negative since
we are trying to find a minimum and therefore descending.


\subsubsection{Overall computation}

Now that we know what we have to compute during each iteration,
let's step back and see how the entire process works.  In the
following algorithm we will assume we find the projection
by using Cholesky factorization to solve the linear system,
as described in Section~\ref{section:projection}.
The linear equation could also be solved using other techniques,
such as iterative solvers.

\begin{enumerate}
\item Find a starting feasible interior point $x_{0}$
\item Symbolically factor $AA^{\tt{T}}$.  Since the non-zero
structure of $BB^{\tt{T}}$ on each step is going to
be the same as $AA^{\tt{T}}$, we can precompute the 
order of the operations in the Cholesky factorization, along
with its non-zero structure.  Since we are only doing this once, we
can potentially spend some time trying to minimize the ``fill'' (the number
of additional non-zero's that appear in the Cholesky factorization that
did not appear in $AA^{\tt{T}}$).
\item Repeat until done 
\begin{itemize}
\item $B = AD_{i}$, where $D_{i}$ is a diagonal matrix with the $x_{i}$ along
      the diagonal.
\item Solve for $w$ in $BB^{\tt{T}}w = BDc$ by multiplying out $BDc$ and
      factoring $BB^{\tt{T}}$ with help from the precomputed symbolic 
      factorization.
\item $z' = - (Dc - B^{\tt{T}}w)$
\item $z = D_{i} z'$  (substituting back to the original space)
\item $x_{i+1} = x_{i} + \alpha z$
\end{itemize}
\end{enumerate}

As mentioned,
we should pick $\alpha$ so that the step is taken some fraction
of the way to the nearest boundary.  See
Figure~\ref{fig:alpha}.  In practice $\alpha$ is selected
so that it goes .96 of the way to the boundary.

\begin{figure}
\centerline{\PSbox{../lecture09/alpha.eps}{1.2in}{2in}}
\caption{Picking an alpha}
\label{fig:alpha}
\end{figure}

\subsection{Potential reduction}
\index{Interior Point: Potential Reduction}

In the potential reduction method, we reword the minimization problem
to be
\begin{quote}
\begin{tabular}{rl}
minimize: & $q\ln (cx -by) - \sum_{j=1}^{n} \ln (x_{j})$ \\
s.t. & $Ax = b$, $x \geq 0$ \\
& $ yA + s = 0$, $s \geq 0 $\\
\end{tabular}
\end{quote}
As we can see, this method makes use of the problem's dual.  The
minimization criteria is partly based on the duality gap ($cx - by$).
Also, the $s$ variables are the slack variables.

The first quantity in the minimization ($q\ln (cx -by)$) represents the
optimizing term.  It will go to $-\infty$ at the optimal solution.
Therefore, this term will dominate near the optimal solution.

$\sum_{j=1}^{n} \ln (x_{j})$ is the analytical center term.

In order to optimize the solution at each iteration, we can either do
a constrained Newton step, or a transformation with a steepest descent.  This
second method is what Karmarkar used in his version of the potential
reduction method.

There are several advantages of this method over the affine scaling
method:
\begin{itemize}
\item There are guaranteed bounds on the running time.
\item The duality gap used in the minimization allows the user to
specify an error tolerance for the solution. 
\item It is more efficient in practice.
\end{itemize}

\subsection{Central trajectory}
\index{Interior Point: Central Trajectory}

The central trajectory is the method of choice these days.  It falls under
the category of Log Barrier methods, which date back to the 50s where
they were used to solve nonlinear problems.

To perform central trajectory, at each iteration step $k$, we do the
following:
\begin{enumerate}
\item 

\begin{tabular}{rl}
minimize: & $cx - \mu_{k}\sum_{j=1}^{n} \ln (x_{j})$ \\
s.t.& $Ax = b, x > 0$ \\
\end{tabular}

(approximate using a constrained Newton step)
\item Select $\mu_{k+1} \leq \mu_{k}$
\end{enumerate}

This will terminate when $\mu_{k}$ approaches 0.

It turns out that with the ``appropriate'' schedule of $\mu_{k}$,
the central trajectory approach reduces to the potential reduction
approach.  We can therefore view it as a generalization.

Currently the primal-dual version of the central trajectory approach
using ``higher order'' approximations (rather than Newton's) seems to
be the best interior-point algorithm in practice.

\section{Algorithm summary}

\begin{enumerate}

\item Overall, the current state of the art uses very sophisticated
algorithms to solve linear programming problems in practice.  Even the
linear algebra packages used to perform the computations use
complicated math in order to optimize computation time.

\item Practice matches theory reasonably well.

\item We find interior-point methods dominant when:
\begin{itemize}
\item Large $n$
\item Small Cholesky factors
\item The problem is highly degenerate, meaning the many of the
constraints cross each other.
\end{itemize}

\item Ellipsoid algorithms are not currently practical, since they always
take the worse case time.

\item Large linear programming problems take hours/days to solve.
Parallelism is very important to speed up the optimization.

\item There are many interior-point methods, but they all share the same
intuition.
\end{enumerate}

We looked at a slide that showed the running times of large
models. For some interior point methods, $AA^ty = Ac$ was too dense, so
the programs ran out of memory and couldn't complete their tasks.

Some algorithms use interior point methods to get close to a solution,
then snap to a boundary and use simplex methods to finish off the
answer.

\section{General Observation}

Mathematics has become increasingly important for ``Real World''
algorithms.  Just over the last 10 years, the complexity of math
employed in algorithms has increased dramatically.  

\begin{itemize}

\item Compression
\begin{enumerate}
\item Information theory
\item Markov models
\item Transforms (for lossy compression)
\end{enumerate}

\item Cryptography
\begin{enumerate}
\item Group theory 
\item Number theory
\end{enumerate}

\item Linear Programming
\begin{enumerate}
\item Matrix algebra
\item Optimization
\item Graph theory
\end{enumerate}
\end{itemize}

