\section{Moving beyond the Simplex Method: An Interior Point \mbox{Algorithm}}\label{S:LP_Interior_Method}
\setcounter{progress_ck_counter}{0}\index{interior point algorithm}

\noi The revised simplex algorithm represents a slight improvement on the original method. However, it does not address a major shortcoming of the simplex algorithm, namely that the number of iterations can grow exponentially large with the size of the problem.
In this section, we investigate a radically different approach for solving LPs.

\subsection{The Origin of the Interior Point Algorithm}\label{S:InteriorMethod_Origin}\index{interior point algorithm!origin}

\noi The \LP LP is a small-scale problem having only two decision variables, three constraints, and two sign restrictions. If we start with the initial basic feasible solution at the origin, we arrive
at the optimal solution after only three iterations of the simplex algorithm. Obviously, for larger scale problems, the number of iterations needed to obtain an optimal solution can become quite large. How much so
 was demonstrated in 1972, when Klee and Minty constructed a family of LPs, each having $n$
 decision variables and $n$ constraints, where $n=1,2,3,\ldots$, yet requiring $2^n$ simplex algorithm iterations to solve \cite{Klee_ExponentialComplexity}. In simple terms, the
 number of iterations needed to solve an LP can be exponential in the problem size. However, this is an extreme case, and for most practical applications, far fewer iterations are required.

\newpar
In 1984, Narendra Karmarkar,\index{Karmarkar, N.} an employee of AT\&T, devised a method for obtaining an approximate
solution to any LP, one that achieved any desired degree accuracy in a number of iterations that was no greater than a polynomial
in the problem size
\cite{Karmarkar_InteriorAlgorithm}, \cite{Dantzig_Thapa}. His method differs fundamentally from the simplex algorithm in that it begins with a feasible solution of the LP having the property that every variable value, be it decision, slack, or excess,
is strictly positive. Such a point is said to belong to the \emph{interior} of the feasible region.
Through an iterative process, a sequence of feasible solutions is constructed, which is contained
in the feasible region's interior and which converges to the optimal solution.

\subsection{The Projected Gradient}\label{S:InteriorMethodProjectedGradient}

\noi Karmarkar's original method is beyond the scope of this text, so we instead focus on a simple variant
known as the \emph{affine scaling algorithm}\index{affine scaling algorithm}. Throughout the discussion, we consider a maximization problem, but to avoid the need for notation that distinguishes among variable types (e.g., decision,
slack, excess), we
express the LP in a form having only equality-type constraints. Namely, we write the
LP as
\begin{align}\label{E:Interior_Method_StandardLPForm}
 \mbox{maximize } z &=\vecc \cdot \vecx  \\
	\mbox{ subject to} \notag\\
	A \vecx &= \vecb \notag \\
	\vecx &\geq \zerovec, \notag
	\end{align}
where $\vecx$ belongs to $\R^n$, $\vecc$ is a row vector in $\R^n$, $\vecb$ belongs to $\R^m$,
and $A$ is an $m$-by-$n$ matrix. For the \LP LP, $\vecx$ belongs to $\R^5$,
\begin{equation*}
\vecc=\begin{bmatrix} 4 & 3 & 0 & 0 & 0 \end{bmatrix} \mbox{,  }
 \vecb = \begin{bmatrix} 8 \\ 28 \\ 32 \end{bmatrix} \mbox{, and }
 A=\begin{bmatrix}
1 & 0 & 1 & 0 & 0 \\
2 & 2 & 0 & 1 & 0 \\
3 & 2 & 0 & 0 & 1 \end{bmatrix}.
\end{equation*}
In this notation, $x_3$, $x_4$, and $x_5$ correspond to the slack variables from the original
\LP LP, and the matrix $A$ denotes the coefficient matrix from the original \LP LP, augmented with
the identity matrix, $I_3$.

\newpar
To begin the algorithm, we start with a point, $\vecx_0$ that belongs to the interior of the LP's feasible region. To say that $\vecx_0$ is feasible means that $A \vecx_0=\vecb$. That
it belongs to the interior means each component of $\vecx_0$ is strictly
positive, as opposed to merely nonnegative. For example, $\ds \vecx_0= \begin{bmatrix} 3 \\ 4 \\ 5 \\ 14 \\15 \end{bmatrix}$ satisfies both these conditions for the \LP LP. Finally, we say
that $\vecx_0$ belongs to the boundary of the feasible region if at least
one of its components is zero.

\newpar
We now start at $\vecx_0$ and seek a new point $\vecx_1=\vecx_0 + \triangle \vecx$, whose
objective, $\vecc \vecx_1$, is larger than at $\vecx_0$. An initial choice for $\triangle \vecx$
is $\vecc^t$, the transpose vector of $\vecc$. That this choice makes sense can be seen if we
observe that $\nabla (\vecc \cdot \vecx) = \vecc^t$, the gradient of $z$, which ``points'' in the
direction of greatest increase in $z$. (We discuss the gradient in much greater detail in Part II of
the text.) However, this choice does not work if we use the feasibility of $\vecx_0$, along
with the desired feasibility of $\vecx_1$. Namely,
\begin{align*}
\vecb &= A\vecx_1 \\
&= A\left(\vecx_0 + \vecc^t \right) \\
&= A\vecx_0 + A \vecc^t \\
&= \vecb + A\vecc^t.
\end{align*}

\newpar
This result implies that $\Delta \vecx=\vecc^t$ must be a solution of the homogeneous matrix equation,
$A\vecx=\zerovec$, or, in other words, $\Delta \vecx=\vecc^t$ belongs to the null space\index{null space} of the
matrix $A$ from (\ref{E:Interior_Method_StandardLPForm}). Such
a requirement does not hold for the case of the \LP LP, nor for any LP
in which the column vectors of $A$ form a linearly independent set.

\newpar
To circumvent this problem, we must therefore choose $\Delta \vecx$ so that it ``points''
as much as possible in the direction of $\vecc^t$, yet belongs to the null space of $A$. One means
of accomplishing this task is to project $\vecc^t$ onto the null space of $A$ by means of a matrix
transformation.

\begin{definition}\label{Th:ProjectedGradient}\index{projected gradient}
Suppose an LP in $m$ constraints and $n$ variables is written in the form
(\ref{E:Interior_Method_StandardLPForm}), where $A$ is an $m$-by-$n$ matrix.
Then the \emph{projected gradient} associated with the LP is defined
as the quantity
\begin{equation}\label{E:Interior_Method_ProjectedGradient}
\vecc_p = P\vecc^t,
\end{equation}
where $P=\left( I_n - A^t\left(AA^t\right)^{-1}A\right)$.
\end{definition}

\newpar
In this definition, we assume that the $m$-by-$m$ matrix $AA^t$ is invertible. This is true
for the \LP LP and any matrix $A$ having linearly independent row vectors.

\begin{progress_ck}
For the \LP LP, calculate the projection matrix $P$\index{projection matrix}. Then use your result
to verify that the projected gradient is given by
\begin{equation*}
\vecc_p=\begin{bmatrix}  \frac{6}{35} \\ \frac{1}{7} \\ -\frac{6}{35} \\ -\frac{22}{25} \\ -\frac{4}{5}
\end{bmatrix}.
\end{equation*}
\end{progress_ck}

\newpar
The projected gradient, $\vecc_p$, and corresponding projection matrix $P$ satisfy three
important properties, whose verifications we leave as exercises:
\be
\item The vector $\vecc_p$ belongs to the null space of $A$. That is, $A\vecc_p=\zerovec$.
\item The matrix $P$ is symmetric, meaning $P^t=P$.
\item The matrix $P$ satisfies $P^2\vecc^t=P\vecc^t$.
\ee
The following theorem is a consequence of these facts.

\begin{theorem}\label{LP_Interior_Method_ImprovingObjective}
For $\alpha >0$, the objective evaluated at $\vecx_1=\vecx_0+\alpha \vecc_p$ is at least as large
as that evaluated at $\vecx_0$. In other words, $\vecc \vecx_1 \geq \vecc \vecx_0$.
\end{theorem}
\begin{proof}
Using the preceding properties, along with properties of the transpose, we
have
\begin{align*}
\vecc \cdot \vecx_1- \vecc \cdot \vecx_0
&= \vecc \left(\vecx_1-\vecx_0 \right) \\
&= \vecc \cdot \alpha \vecc_p \\
&= \alpha \vecc \left(P\vecc^t\right) \\
&= \alpha \vecc \left(P^2 \vecc^t \right) \\
&= \alpha \left(\vecc P^t \right)\left(P \vecc^t\right) \\
&= \alpha \left(P \vecc^t \right)^t \left(P \vecc^t\right) \\
&= \alpha \left\| \vecc_p \right\|^2 \\
&\geq 0.
\end{align*}
\end{proof}
\subsection{Affine Scaling}\label{S:InteriorMethod_AffineScaling}\index{affine scaling}

\noi Theorem \ref{LP_Interior_Method_ImprovingObjective} demonstrates
that for $\alpha >0$, $\vecx_1=\vecx_0+\alpha \vecc_p$ has an objective value at least as large as that of
$\vecx_0$ and that this value increases as $\alpha$ increases. Since $\vecc_p$ belongs to the null space of $A$,
 $A\vecx_1=A \vecx_0+\alpha A\vecc_p = A \vecx_0 \vecb$. However, this fact alone does not guarantee
 that $\vecx_1$ is feasible. For if $\alpha$ increases by too much, there exists the possibility
 that a component of $\vecx_1$ becomes negative, in which case $\vecx_1$ does not
 satisfy the sign restrictions. The key then is to allow $\alpha$ to increase, but not by too
 much.

\newpar
The \LP LP illustrates the delicacy of this situation. Using
$\ds \vecx_0= \begin{bmatrix} 3 \\ 4 \\ 5 \\ 14 \\15 \end{bmatrix}$ and the result
from the previous Waypoint, we have
\begin{align}\label{E:LP_Interior_Method_Fuelpro_alphachoice}
\vecx_1 &= \vecx_0 + \alpha \vecc_p \\
&= \begin{bmatrix} 3 \\ 4 \\ 5 \\ 14 \\15 \end{bmatrix} + \alpha
\begin{bmatrix}  \frac{6}{35} \\ \frac{1}{7} \\ -\frac{6}{35} \\ -\frac{22}{25} \\ -\frac{4}{5}
\end{bmatrix} \notag \\
&= \begin{bmatrix}
3 + \frac{6}{35}\alpha  \\
4 + \frac{1}{7}\alpha\\
5  -\frac{6}{35}\alpha \\
14  -\frac{22}{25}\alpha \\
15 -\frac{4}{5}\alpha
\end{bmatrix}. \notag
\end{align}

\noi
The last component of this vector dictates that $\alpha$ can increase by
no more than $18 \frac{3}{4}$ before $\vecx_1$ violates the sign restrictions.

\newpar
Setting $\ds \alpha= 18 \frac{3}{4}$ results in
\begin{equation}\label{E:LP_Interior_Method_Fuelpro_alphabinding}
\vecx_1 = \begin{bmatrix}
\frac{87}{14} \\
\frac{187}{28} \\
\frac{25}{14} \\
\frac{31}{14} \\
0
\end{bmatrix},
\end{equation}
with corresponding objective value $z=\vecc \cdot \vecx_1=\frac{1257}{28} \approx 44.9$.

\newpar
This objective value is much larger than $z=\vecc \cdot \vecx_0=24$ and much
closer to the known optimal value of $z=46$. Moreover, the last component
of $\vecx_1$ is zero, implying that $\vecx_1$ corresponds to a point
on the boundary of the original LP's feasible region where
the third constraint, $3x_1+2x_2 \leq 32$, is binding. However,
the values of the original LP's decision variables, $\ds x_1 =\frac{87}{14} \approx 6.21$,
and $\ds x_2 = \frac{187}{28} \approx 6.67$, are far from close
to the known optimal solution, $x_1=4$ and $x_2=10$.

\newpar
Repeating this entire process by computing a new value $\vecx_2$ from $\vecx_1$ in the same manner
we computed $\vecx_1$ from $\vecx_0$ will not improve the quality of the estimate. This follows from
observing that the last component of $\vecx_2=\vecx_1 + \alpha \vecc_p$ is
$\ds -\frac{4}{5}\alpha$. For $\vecx_2$ to remain feasible, $\alpha$ must equal zero, whence
$\vecx_2=\vecx_1$.

\newpar
This result demonstrates a fundamental problem we must overcome, which can be summarized in the following terms. The quantity $\alpha$
can be viewed as a ``step size,'' which measures what fraction of the projected gradient,
$\vecc_p$, we use to obtain the next approximation, $\vecx_{k+1}$, from the current one, $\vecx_k$.
On the one hand, we desire a large step size so that the difference in objective values,
$\ds \vecc \vecx_{k+1}-\vecc \vecx_k = \alpha \vecc \cdot \vecc_p$, is significant.
On the other hand, if the step size, $\alpha$, is too large, a component of
$\vecx_{k+1}=\vecx_k + \alpha \vecc_p$ becomes negative so that $\vecx_{k+1}$ violates
the sign restrictions.

\newpar
The key to overcoming this problem is to perform a change of variables at each iteration by means of a simple
matrix transformation. To compute
$\vecx_{k+1}$, we first use the previous iterates, $\vecx_0,\vecx_1,\ldots,\vecx_k$,
to create a completely new LP. In this new LP, the initial value, $\tilde{\vecx_0}$, which
corresponds to $\vecx_k$ in the original LP,
lies far enough from the boundary
of the feasible region so that
if we set it to the initial value in new LP, we are permitted to select a large step size. After using this large step size to calculate a better approximate solution to the new LP, we ``change back'' variables, so to speak, to determine the corresponding value
of $\vecx_{k+1}$. This overall strategy is illustrated in Figure \ref{Fi:InteriorPointAlgorithmCommDiag}.

\scalebox{.95}{\begin{minipage}{\linewidth}
\begin{equation*}
\begin{CD}
\mbox{New LP: max } z=\vecc \cdot D\tilde{\vecx} \\
\mbox{such that }AD\tilde{\vecx}=\vecb \mbox{, }
\tilde{\vecx} \geq \zerovec \\
\mbox{Set } \tilde{\vecx}_0=D^{-1}\vecx_k
 @>>>
 \begin{matrix}
 \mbox{Compute }  \tilde{\vecx}_1 \mbox{ via} \\
 \mbox{projected gradient of } AD
\end{matrix}   \\
@AA A     @VVV    \\
@AA  \tilde{\vecx}=D^{-1}\vecx      A     @V\vecx=D\tilde{\vecx} VV    \\
@AA\text{Centering Transformation}     A     @V\text{Transform back to original LP}VV   \\
@AA A     @VVV    \\
\mbox{Original LP: max } z=\vecc \cdot \vecx \\
\mbox{such that }A\vecx=\vecb \mbox{, }
\vecx \geq \zerovec \\
\mbox{starting with } \vecx_k
 @<\mbox{Repeat the process}<<
\vecx_{k+1}=D\tilde{\vecx}_1
\end{CD}
\end{equation*}

\end{minipage}}

\vspace{-.2in}
\begin{figure}[h]
\caption{Diagram illustrating change of variables process.}\label{Fi:InteriorPointAlgorithmCommDiag}
\end{figure}

\newpage
\noindent The matrix $D$ that we elect to use has the effect of ``centering''
$\tilde{\vecx_k}$ in the feasible region of the new LP. It is an $n$-by-$n$ matrix defined by
\begin{equation}\label{E:LP_Interior_Method_CenteringMatrix}
D= \begin{bmatrix}
x_{k,1} & 0 & \ldots & 0 \\
0 & x_{k,2}& \ldots  & 0  \\
\vdots & & \ddots & \vdots  \\
0 & 0 & \ldots & x_{k,n}
\end{bmatrix},
\end{equation}

where $x_{k,i}, 1 \leq i \leq n$ denotes the $i$th component of $\vecx_k$.

\newpar
Here are two important facts regarding the matrix $D$.
\bi
\item Since every component of $\vecx_k$ is strictly positive, $D$ is invertible and
\begin{equation*}
D^{-1} =
\begin{bmatrix}
x_{k,1}^{-1} & 0 & \ldots & 0 \\
0 & x_{k,2}^{-1}& \ldots  & 0  \\
\vdots & & \ddots & \vdots  \\
0 & 0 & \ldots & x_{k,n}^{-1}
\end{bmatrix}
\end{equation*}.
\item In the new LP,
each component of $\tilde{\vecx_0}=D^{-1}\vecx_k$ equals one. This means that each component of $\tilde{\vecx_0}$ is one unit from the boundary of the interior of the
    new LP. To emphasize this crucial fact, we use the notation, $\onevec_n$ to denote the
    vector $\tilde{\vecx_0}$.
\ei

\newpar
\subsection{Summary of the Method}\label{S:InteriorMethod_Outline}\index{interior point algorithm!summary of steps}

\noi Here then is an overview of our interior point algorithm.
\be
\item Suppose $\vecx_k$, where $k \geq 0$, lies in the interior of the feasible region of
the original LP, meaning $A\vecx_k=\vecb$ and each component of $\vecx_k$ is strictly positive.
\item We construct a diagonal matrix, $D$, whose diagonal values are the components of $\vecx_k$.
\item In the new LP, the objective coefficient is given by $\tilde{\vecc}=\vecc D$ and the
coefficient matrix by $\tilde{A}=AD$. We compute the projected gradient of $\tilde{\vecc}$ for the new LP. That
is, we project $\tilde{\vecc}$ onto the null space of $\tilde{A}$ and set
\begin{equation}\label{E:LP_Interior_Method_NewProjectedGradient}
\tilde{\vecc_p} = P\tilde{\vecc}^t,
\end{equation}
where
\begin{equation*}
P=\left( I_n - \tilde{A}^t\left(\tilde{A}\tilde{A}^t\right)^{-1}\tilde{A}\right).
\end{equation*}

\item Since $\tilde{\vecx_0}= \vece$, where $\vece$ is the vector in $\R^n$ all of
whose entries are one, $\alpha$ can be at most the reciprocal of the absolute value of the most negative
entry of $\tilde{\vecc_p}$ before an entry of
$\tilde{\vecx_0}+\alpha \tilde{\vecc_p}$ becomes negative. We therefore
set $\alpha$ to equal this quantity.

 \item Define
\begin{equation*}
\tilde{\vecx_1} = \vece + \lambda \alpha  \tilde{\vecc_p},
\end{equation*}
where $\lambda$ is a fixed quantity slightly less than one. The value of $\lambda\alpha$
controls the stepsize used in the new LP.

\item Let $\vecx_{k+1}=D\tilde{\vecx_1}$, the $(k+1)$st iterate in the original LP.
\item Return to Step 1 and repeat the process.
\ee

\subsection{Application of the Method to the \LP LP}\label{S:InteriorMethod_Fuelpro}

\noi We now demonstrate the details for performing one iteration of the algorithm
for the \LP LP.
\be
\item We start with the initial value
$\ds \vecx_0= \begin{bmatrix} 3 \\ 4 \\ 5 \\ 14 \\15 \end{bmatrix}$.
\item Using $\vecx_0$,
\begin{equation*}
D=\begin{bmatrix}
3 & 0 & 0 & 0 & 0\\
0 & 4 & 0 & 0 & 0\\
0 & 0 & 5 & 0 & 0\\
0 & 0 & 0 & 14 & 0\\
0 & 0 & 0 & 0 & 15\\
\end{bmatrix}
\end{equation*}.

\item In the new LP, $\ds \tilde{\vecc}=\vecc D = \begin{bmatrix} 12 & 12 & 0 & 0 & 0 \end{bmatrix}$
and \newline
$\tilde{A}= AD= \begin{bmatrix} 3 & 0 & 5 & 0 & 0 \\
6 & 8 & 0 & 14 & 0 \\
9 & 8 & 0 & 0 & 15 \end{bmatrix}$.
These quantities are then used to construct the projection matrix, $P$, and projected gradient,
$\tilde{\vecc_p}$. The
second of these is given by
$\ds \tilde{\vecc_p} \approx \begin{bmatrix} 4.57 \\ 5.85 \\ -2.74 \\ -5.30 \\ -5.86 \end{bmatrix}$.
\item The last entry of $\tilde{\vecc_p}$ indicates that $\alpha=1/5.86 \approx .1706$, the reciprocal of the
absolute value of most negative component of $\tilde{\vecc_p}$. Using this value along
with $\lambda=.75$, we compute $\tilde{\vecx_1}$:
\begin{align*}
\tilde{\vecx_1} &= \tilde{\vecx_0} + \lambda \alpha  \tilde{\vecc_p} \\
& \approx \begin{bmatrix} 1.58 \\ 1.75 \\ .65 \\ .32 \\ .25 \end{bmatrix}
\end{align*}
\item In the original LP,
\begin{equation*}
\vecx_{1} = D \tilde{\vecx_1} \approx \begin{bmatrix} 4.75 \\ 6.99 \\ 3.25 \\ 4.50 \\ 3.75 \end{bmatrix},
\end{equation*}
which has corresponding objective value, $\vecc \cdot \vecx_1 \approx 40$.
\ee

\newpar
\subsection{A Maple Implementation of the Interior Point \mbox{Algorithm}}\label{SS:InteriorMethod_Maple}\index{interior point algorithm! Maple implementation}

\noi The affine scaling method described in this section can be implemented quite easily using Maple. What follows is a worksheet, \textbf{Interior Point Algorithm.mw}, that computes the four iterations.

\begin{maplelist}
\item \mf{restart;with(LinearAlgebra):}
\item \mf{c:=Vector[row]([4,3,0,0,0]); \newline \# Create row vector with objective coefficients.}
\mo{c:=[4,3,0,0,0]}
\item \mf{A:=Matrix(3,5,[1,0,1,0,0, 2,2,0,1,0, 3,2,0,0,1]); \newline \# Matrix of constraint coefficients.}
\mo{A:=\begin{bmatrix} 1 & 0 & 1 & 0 & 0 \\ 2 & 2 & 0 & 1 & 0 \\
3 & 2 & 0 & 0 & 1 \end{bmatrix}}
\item \mf{b:=<8,28,32>; \newline \# Constraint bounds.}
\mo{b:=\begin{bmatrix} 8 \\ 28 \\ 32 \end{bmatrix} }
\item \mf{N:=4:lambda:=.75: \newline \# Set number of iterations, N, and parameter lambda.}
\item \mf{x:=array(0..N): \newline \# Array of iteration values.}
\item \mf{x[0]:=<3,4,5,14,15>;\newline \# Set initial value.}
\mo{x_0:=\begin{bmatrix} 3 \\ 4 \\ 15 \\ 14 \\ 15 \end{bmatrix}}
\item \mf{for i from 0 to (N-1) do \newline
d:=DiagonalMatrix(convert(x[i],list)):
\newline \# Determine transformation matrix. \newline
$\sim$ A:= A.d: $\sim$ c:=c.d:
\newline \# Calculate coefficient matrix and objective coefficients
\newline for new LP.
P:=IdentityMatrix(5)- \newline Transpose($\sim$ A).MatrixInverse($\sim$
A.Transpose($\sim$ A)).$\sim$ A:
\newline \# The projection matrix. \newline
cp:=P.Transpose($\sim$ c):
\newline \# Determine projected gradient. \newline
alpha:=(abs(min(seq(cp[j],j=1..5))))\^{}(-1):
\newline \# Find alpha. \newline
x[i+1]:=d.(<1,1,1,1,1>+lambda*alpha*cp):\newline
\# Determine next iterate. \newline
end do:}
\item \mf{for i from 0 to N do print(x[i],Transpose(c).x[i]):end do:
\newline \# Print sequence of iterates, together with corresponding 
\newline objective values.}

\end{maplelist}
For the sake of brevity, we do not print here the output that appears at the end of this worksheet. Instead, we summarize these results in Table \ref{Ta:LP_Interior_Method_IterationTable} but list only the entries of each $\vecx_i$ that
correspond to the two decision variables, $x_1$ and $x_2$, in the original \LP LP.

\begin{table}[h]
\caption{Results of applying the interior point algorithm to the \LP LP}
\label{Ta:LP_Interior_Method_IterationTable}
\begin{tabular}{|c|c|c|c|}
\hline
$n$ &  $(x_{n,1},x_{n,2}) $  & z \\
\hline
0 & $(3,4)$ &    24 \\ \hline
1 & $(4.75,7) $ &  40   \\ \hline
2 & $(5.1, 7.95) $ &   44.07  \\ \hline
3 & $(4.66, 8.9) $ &   45.32  \\ \hline
4 & $(4.1, 9.79)$ &   45.76  \\ \hline
\end{tabular}
\end{table}	

\newpar
The results in Table \ref{Ta:LP_Interior_Method_IterationTable} dictate that
convergence to the optimal solution takes place very quickly. A graphical illustration
of this phenomenon is depicted in Figure \ref{Fi:f03x06x01}.

\begin{figure}[h]
\centering
\psfrag{x0}[Bc][Bc]{$\vecx_0$}
\psfrag{x1}[Bc][Bc]{$\vecx_1$}
\psfrag{x2}[Bc][Bc]{$\vecx_2$}
\psfrag{x3}[Bc][Bc]{$\vecx_3$}
\psfrag{x4}[Bc][Bc]{$\vecx_4$}
\includegraphics[width=.7\linewidth]{images/f03x06x01.eps}
\caption{Interior algorithm iterates for \LP LP.}\label{Fi:f03x06x01}
\end{figure}

\newpar
Since 1984, interior point methods have continually evolved and improved in
their efficiency. Interestingly, after a great deal of controversy centered on the issuance of patents for what could be viewed as newly developed mathematical tools, AT\& T applied for, and received a patent for Karmarkar's method in 1988. This patent expired in 2006.

\section*{Exercises Section \protect\ref{S:LP_Interior_Method}}

\addcontentsline{toc}{subsection}{Exercises}

\input{Exercises/Exercises_LP_Interior_Method}		
\endinput
	
	