% $Header: /home/rdv/quantum/genetic-paper/RCS/eaqc-paper.tex,v 1.12 2005/05/29 21:34:00 rdv Exp rdv $
% now how do I set the paper size?
% \documentclass[twocolumn]{article}

\documentclass[onecolumn,doublespace]{article}
\usepackage{times}
\usepackage{graphics}
\usepackage{lscape}
% \usepackage{subfigure}
\usepackage{caption}
\usepackage{subcaption}
\usepackage{graphicx}
\usepackage{color}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{kbordermatrix,blkarray}
\pagestyle{plain}

\bibliographystyle{unsrt}   %>>>> makes bibtex use spiebib.bst

% \input{psfig}
\input{epsf}
% \epsfverbosetrue

%\def\comment#1{$\Rightarrow${\small\em #1}}
\def\comment#1{}
\newenvironment{chapterquote}
  {\begin{quotation}
      \sf}
  {\end{quotation}}
\def\cq#1#2{\begin{chapterquote}#1\flushright\textbf{#2}\end{chapterquote}}


\newcommand{\icol}[1]{% inline column vector
  \left(\begin{smallmatrix}#1\end{smallmatrix}\right)%
}

\newcommand{\irow}[1]{% inline row vector
  \begin{smallmatrix}(#1)\end{smallmatrix}%
}

\newcommand{\tvec}[1]{%
  \vec{#1}^{\thinspace\intercal}%
}

\begin{document}
\title{Optimization Theory (DS2) Lecture \#4\\
The Heart of the Simplex Algorithm
}
\maketitle
\begin{abstract}
  Today we are going to look at the simplex algorithm, as developed by
  George Dantzig (1947).  We will cover Secs. 2.4 and 2.5 of the
  textbook.  \emph{In my opinion, this is the hardest week of the
    semester, stick with me here!}
\end{abstract}

\emph{Text and examples adapted from \emph{A Gentle Introduction to
    Optimization}.}

\section{Review}

There were a lot of blank looks last time as we discussed the
geometric nature of the constraints and the objective function, and we
didn't even get to discuss \emph{standard equality form} (SEF) and the
core idea of the \emph{simplex operation}.  So let's begin here with
some geometry, then go back and pick up what we were supposed to cover
in Sec. 4 and 5 of last time, before we start into the grungy parts of
simplex.

\subsection{The Geometric Nature of Constraints and Objective
  Functions}

Or, \emph{A Simple(!) Example of Prepping a Linear Program on a 2-D
  Simplex (i.e., a Triangle)}

By now, you've heard me say that a \emph{simplex} is an
$n$-dimensional triangle, but that the \emph{simplex algorithm}
applies to \emph{polytopes} ($n$-dimensional polygons).

Let's look at a problem formulation for a 2-D simplex, that is, a
triangle.  Examples in this section are in Octave, rather than in R.

Consider the following LP, not yet in SEF (defined in the next
section):

\begin{equation}
  \operatorname{max} z(\vec{x}) = (\frac{19}{14},-\frac{1}{28})(x_1,x_2)^\intercal
\end{equation}
subject to
\begin{align}
  \left(
    \begin{array}{cc}
      1 & 0 \\
      0 & 1 \\
      2 & 3
    \end{array}
  \right)
  \vec{x}
  \begin{array}{c}
    \ge \\
    \ge \\
    \le
  \end{array}
  \left(
    \begin{array}{c}
      2 \\
      4 \\
      28
    \end{array}
  \right) \\
  x_1, x_2 \ge 0.
\end{align}

It should be obvious that the first constraint is $x_1 \ge 2$, and the
second constraint is $x_2 \ge 4$.  The third constraint, $2x_1 + 3x_2
\le 28$, is a diagonal line.  These constraints, one at a time, then
the whole set, are illustrated in Fig.~\ref{fig:constraints}.  Our
objective function is shown in Fig.~\ref{fig:triangle-objective}, with
the code for all of these figures in the next subsection.

\begin{figure}[h]
  \begin{subfigure}{0.5\textwidth}
    \includegraphics[width=\textwidth]{constraint1}
    \caption{$x_1 \ge 2$}
    \label{fig:constraint1}
  \end{subfigure}
  \begin{subfigure}{0.5\textwidth}
    \includegraphics[width=\textwidth]{constraint2}
    \caption{$x_2 \ge 4$}
    \label{fig:constraint2}
  \end{subfigure}
  \begin{subfigure}{0.5\textwidth}
    \includegraphics[width=\textwidth]{constraint3}
    \caption{$2x_1 + 3x_2 \le 28$}
    \label{fig:constraint3}
  \end{subfigure}
  \begin{subfigure}{0.5\textwidth}
    \includegraphics[width=\textwidth]{constraints-all}
    \caption{All three constraints.}
    \label{fig:constraints-all}
  \end{subfigure}
  \caption{The constraints in our simplex example.  The three
    constraints collectively define a triangle of acceptable solutions.}
  \label{fig:constraints}
\end{figure}

\begin{figure}
\includegraphics[width=13cm]{triangle-objective.pdf}
\caption{The objective function $z(\vec{x}) =
  \frac{19}{14}x_1-\frac{1}{28}x_2$ in our 2-D simplex, bounded by the
  constraints above.  For linear programs (problems), the optimal
  point always lies at a vertex (corner).}
\label{fig:triangle-objective}
\end{figure}
Our first important task is to put the problem into standard equality
form (SEF).  Since our constraints are inequalities, we have to add
slack variables, $x_3$, $x_4$ and $x_5$.  (Slack variables and free
variables are covered more rigorously in Sec.~\ref{sec:sef}.)

\begin{equation}
  \operatorname{max} z(\vec{x}) = (\frac{19}{14},-\frac{1}{28})(x_1,x_2)^\intercal
\end{equation}
subject to
\begin{align}
  \left(
    \begin{array}{ccccc}
      1 & 0 & -1 & 0 & 0 \\
      0 & 1 & 0 & -1 & 0 \\
      2 & 3 & 0 & 0 & 1
    \end{array}
  \right)
  \vec{x} & =
  \left(
    \begin{array}{c}
      2 \\
      4 \\
      28
    \end{array}
  \right) \\
  x_1, x_2, x_3, x_4, x_5 & \ge 0.
\end{align}

Now our program is in SEF, nice and pretty and waiting to be solved!

\subsection{Interlude in Octave}

We can draw each of those constraints using the following Octave code:

\begin{verbatim}
>> c1 = [0 0; 0 10; 2 10; 2 0];
>> c1v = [0 0; 0 10; 2 10; 2 0];
>> c2v = [0 0; 0 4; 10 4; 10 0];
>> fs = [ 1 2 3 4];
>> c3v = [0 28/3; 0 10; 10 10; 10 8/3];
>> fs3 = [1 2 3 4];
>> figure
>> xlim([0 10])
>> ylim([0 10])
>> axis('equal')
>> patch("Faces", fs, "Vertices", c1v, 'FaceColor', 'blue');
>> figure
>> xlim([0 10])
>> ylim([0 10])
>> axis('equal')
>> patch("Faces", fs, "Vertices", c2v, 'FaceColor', 'yellow');
>> figure
>> xlim([0 10])
>> ylim([0 10])
>> axis('equal')
>> patch("Faces", fs3, "Vertices", c3v, 'FaceColor', 'red');
>> figure
>> xlim([0 10]); ylim([0 10]); axis('equal');
>> patch("Faces", fs, "Vertices", c1v, 'FaceColor', 'blue');
>> patch("Faces", fs, "Vertices", c2v, 'FaceColor', 'yellow');
>> patch("Faces", fs3, "Vertices", c3v, 'FaceColor', 'red');
\end{verbatim}

And plot the constrained objective function using:

\begin{verbatim}
>> v = [2 4; 2 8; 8 4]
v =

   2   4
   2   8
   8   4

>> f = [1 2 3]
f =

   1   2   3
>> col = [2.57; 1; 5]
col =

   2.5700
   1.0000
   5.0000

>> figure
>> xlim([0 10])
>> ylim([0 10])
>> axis('square')
>> axis('equal')
>> patch('Faces',f,'Vertices',v,'FaceVertexCData',col,'FaceColor','interp');
>> colorbar
\end{verbatim}

\subsection{Standard Equality Form (SEF)}
\label{sec:sef}

\emph{Adapted from Sec. 4 of last time's text.}

An LP is in \emph{standard equality form} if:
\begin{enumerate}
\item it is a maximization problem;
\item other than the non-negativity constraints, all constraints are
  \emph{equalities}; and
\item every variable has a non-negativity constraint.
\end{enumerate}
That is, it can be written in the form
\begin{equation}
  \operatorname{max}\{z(\vec{x}) = \tvec{c}\vec{x} : A\vec{x} =
  \vec{b}, \vec{x} \ge \vec{0}\}.
\end{equation}

In this formulation, $\vec{c}$ is the vector of coefficients of each
of our variables in the objective function.  It will have two entries
for each free variable we adapted and will have zeroes for each slack
variable we added.  $A$ and $\vec{b}$ have one line for each
constraint in our problem.

Consider the LP
\begin{equation}
  \operatorname{max}(1,-2,4)(x_1,x_2,x_3)^\intercal
\end{equation}
subject to
\begin{align}
  \left(
    \begin{array}{ccc}
      1 & 5 & 3 \\
      2 & -1 & 2 \\
      1 & 2 & -1 
    \end{array}
  \right)
\left(
  \begin{array}{c}
    x_1 \\
    x_2 \\
    x_3
  \end{array}
  \right)
  \begin{array}{c}
    \ge \\
    \le \\
    =
  \end{array}
\left(
  \begin{array}{c}
    5 \\
    4 \\
    2
  \end{array}
  \right)\\
x_1, x_2 \ge 0.
\end{align}

$x_3$ is called a \emph{free variable} since it doesn't have a
non-negativity constraint.  But a lot of our techniques and especially
proofs depend on that.  So, we introduce two new variables, $x_3^+$
and $x_3^-$, set $x_3 = x_3^+ - x_3^-$, and add constraints $x_3^+,
x_3^- \ge 0$.  After some algebra, you get
\begin{equation}
  \operatorname{max}(1,-2,4,-4)(x_1,x_2,x_3^+,x_3^-)^\intercal
\end{equation}
subject to
\begin{align}
  \left(
    \begin{array}{cccc}
      1 & 5 & 3 & -3 \\
      2 & -1 & 2 & -2 \\
      1 & 2 & -1 & 1 
    \end{array}
  \right)
\left(
  \begin{array}{c}
    x_1 \\
    x_2 \\
    x_3^+ \\
    x_3^-
  \end{array}
  \right)
  \begin{array}{c}
    \ge \\
    \le \\
    =
  \end{array}
\left(
  \begin{array}{c}
    5 \\
    4 \\
    2
  \end{array}
  \right) \label{eq:sef-pre}\\
x_1, x_2, x_3^+, x_3^- \ge 0.
\end{align}

Getting closer, but we're still not quite there; we don't have
equalities everywhere yet.  So we introduce two \emph{slack variables}
$x_4$ and $x_5$.  Now we get

\begin{equation}
  \operatorname{max}(1,-2,4,-4,0,0)(x_1,x_2,x_3^+,x_3^-,x_4,x_5)^\intercal
\end{equation}
subject to
\begin{align}
  \left(
    \begin{array}{cccccc}
      1 & 5 & 3 & -3 & 0 & -1\\
      2 & -1 & 2 & -2 & 1 & 0 \\
      1 & 2 & -1 & 1 & 0 & 0
    \end{array}
  \right)
\left(
  \begin{array}{c}
    x_1 \\
    x_2 \\
    x_3^+ \\
    x_3^- \\
    x_4 \\
    x_5
  \end{array}
  \right)
    =
\left(
  \begin{array}{c}
    5 \\
    4 \\
    2
  \end{array}
  \right)\\
x_1, x_2, x_3^+, x_3^-, x_4, x_5 \ge 0
\end{align}
and we're finally in SEF!  Note that when we replaced a $\le$
constraint (second line of Eq.~\ref{eq:sef-pre}), we added a $1$ to
the array, and when we replaced a $\ge$ constraint (first line of
Eq.~\ref{eq:sef-pre}), we added $-1$ to the array.  This helps us keep
our non-negativity constraint for all variables.

\subsection{A Simplex Iteration}

\emph{Adapted from Sec. 5 of last time's text.}

The basic idea is to increase one of our variables $x_i$ until it hits
a ``stop'' against one of the constraints.  In the main body of the
simplex algorithm, this will involve running along an edge of the
\emph{polytope} (see below) until we reach a vertex.

Consider the following LP in SEF:

\begin{equation}
  \operatorname{max} z(\vec{x}) = (2,3,0,0,0)(x_1,x_2,x_3,x_4,x_5)^\intercal
  \label{eq:fav-obj-1st-time}
\end{equation}
subject to
\begin{align}
  \left(
    \begin{array}{ccccc}
      1 & 1 & 1 & 0 & 0 \\
      2 & 1 & 0 & 1 & 0 \\
      -1 & 1 & 0 & 0 & 1 
    \end{array}
  \right)
  \vec{x} =
  \left(
    \begin{array}{c}
      6 \\
      10 \\
      4
    \end{array}
  \right) \\
  x_1, x_2, x_3, x_4, x_5 \ge 0.
  \label{eq:fav-nonneg-1st-time}
\end{align}

Given the solution $\vec{x_a} = (0,0,6,10,4)^\intercal$ (easy to see that's
a solution -- look at the right side of the array), $z(\vec{x_a}) =
0$.

Now try increasing $x_1$, choosing $x_1 = t$.  A little algebra gives
\begin{equation}
\left(
  \begin{array}{c}
    x_3 \\
    x_4 \\
    x_5
  \end{array}
  \right) 
= \left(
  \begin{array}{c}
    6 \\
    10 \\
    4
  \end{array}
  \right) 
- t \left(
  \begin{array}{c}
    1 \\
    2 \\
    -1
  \end{array}
  \right) 
\end{equation}
from which we get
\begin{equation}
  t
\left(
  \begin{array}{c}
    1 \\
    2 \\
    -1
  \end{array}
  \right) 
\le \left(
  \begin{array}{c}
    6 \\
    10 \\
    4
  \end{array}
  \right) 
\end{equation}
which gives us the limiting inequality $t \le 10/2 = 5$.  Setting $t = 5$,
$\vec{x'} = (5, 0, 1, 0, 9)^\intercal$ and $z(\vec{x'}) = 10$.

Unfortunately, we can't yet apply the same trick to $x_2$, so that's
the topic for this week.

\section{Bases}

(Adapted from Sec. 2.4.1 in the textbook.)

Before we get to the algorithm itself, we need to discuss the idea of
a \emph{basis} for a set of variables (dimensions), and
\emph{canonical form}.

Recall that the \emph{inverse} of a square matrix $A$ is the matrix
$A^{-1}$ such that $AA^{-1} = I$, where $I$ is the identity matrix.  A
\emph{singular} matrix has no inverse, or equivalently, its
determinant is 0.  A \emph{nonsingular} matrix has an inverse and a
non-zero determinant.  (Question: is the identity matrix singular or
nonsingular?)

$A$ is the $m\times n$ matrix with $m$ linearly independent rows that
comprise our constraints.  We need to partition our columns into two
sets, one that comprises a \emph{basis}, which we will call $B$, and
the rest, which we will call $N$.

Consider

\begin{align}
A = \left(
    \begin{array}{cccccc}
      2 & 1 & 2 & -1 & 0 & 0 \\
      1 & 0 & -1 & 2 & 1 & 0 \\
      3 & 0 & 3 & 1 & 0 & 1
    \end{array}
  \right).
\end{align}
The set of columns $B = \{2, 5, 6\}$ is a basis, which we can see
easily by taking
\begin{align}
A_B = \left(
    \begin{array}{ccc}
      1 & 0 & 0 \\
      0 & 1 & 0 \\
      0 & 0 & 1
    \end{array}
  \right),
\end{align}
which, if you think about it as a three-dimensional space, is the
vector composed of the $x$, $y$ and $z$ unit vectors.

\emph{Question: For the simplex algorithm, does our basis set have to
  be orthonormal?}

The set of columns $B = \{1, 5, 6\}$ is also a basis:
\begin{align}
A_B = \left(
    \begin{array}{ccc}
      2 & 0 & 0 \\
      1 & 1 & 0 \\
      3 & 0 & 1
    \end{array}
  \right).
\end{align}
You can see pretty easily how to turn that into the identity matrix by
substracting $0.5$ times the first row from the second, and $1.5$
times the first row from the third, after which our basis is
orthonormal.

We will need a \emph{basic solution} $\vec{x_a}$ for our problem
$A\vec{x} = \vec{b}$.  $\vec{x_a}$ is a basic solution iff:
\begin{enumerate}
\item $A\vec{x_a} = \vec{b}$, and
\item $x_N = \vec{0}$.
\end{enumerate}
where $x_N$ is the subset of the vector $\vec{x_a}$ of the columns
\emph{not} in the basis $B$.

Every basis has a unique basic solution, which we are going to need.
For  $B = \{1, 5, 6\}$, the basic solution is $\vec{x_a} = (1, 0, 0,
0, 0, -2)^\intercal$.

\section{Canonical Forms}

(Adapted from Sec. 2.4.2 in the textbook.)

Consider an LP (P) in SEF:
\begin{equation}
  \operatorname{max}\{\tvec{c}\vec{x} + z_a: A\vec{x} = \vec{b}, \vec{x}
  \ge \vec{0}\}
\end{equation}
where $z_a$ is a constant.  Let $B$ be a basis of $A$.  (P) is
in \emph{canonical form for $B$} if:
\begin{enumerate}
\item $A_B$ is an identity matrix, and
\item $\vec{c_B} = \vec{0}$.
\end{enumerate}

Let's return to the LP in
Eqs.~\ref{eq:fav-obj-1st-time}-\ref{eq:fav-nonneg-1st-time}
above:


\begin{equation}
  \operatorname{max} z(\vec{x}) = (2,3,0,0,0)(x_1,x_2,x_3,x_4,x_5)^\intercal
  \label{eq:fav-obj-2nd-time}
\end{equation}
subject to
\begin{align}
  \left(
    \begin{array}{ccccc}
      1 & 1 & 1 & 0 & 0 \\
      2 & 1 & 0 & 1 & 0 \\
      -1 & 1 & 0 & 0 & 1 
    \end{array}
  \right)
  \vec{x} =
  \left(
    \begin{array}{c}
      6 \\
      10 \\
      4
    \end{array}
  \right) \\
  x_1, x_2, x_3, x_4, x_5 \ge 0.
  \label{eq:fav-nonneg-2nd-time}
\end{align}

Note that in this case, $z_a = 0$.

Obviously, $B = \{ 3, 4, 5 \}$ is a basis, but that's boring,
it's already in canonical form.  $B = \{ 1, 2, 4 \}$ is also a basis,
so let's rewrite our program in canonical form for this basis:

\begin{equation}
  \operatorname{max} z(\vec{x}) = 17 + (0,0,-5/2,0,-1/2)(x_1,x_2,x_3,x_4,x_5)^\intercal
  \label{eq:fav-basis-124-obj}
\end{equation}
subject to
\begin{align}
  \left(
    \begin{array}{ccccc}
      1 & 0 & 1/2 & 0 & -1/2 \\
      0 & 1 & 1/2 & 0 & 1/2 \\
      0 & 0 & -3/2 & 1 & 1/2 
    \end{array}
  \right)
  \vec{x} =
\left(
  \begin{array}{c}
    1 \\
    5 \\
    3
  \end{array}
  \right) \\
x_1, x_2, x_3, x_4, x_5 \ge 0.
  \label{eq:fav-basis-124-nonneg}
\end{align}

(Notice that now columns 1, 2, 4 form an identity matrix.)  So how did
we get here?  First, we are going to need the \emph{inverse} of a
matrix or two, so let's take a quick break and see how to calculate
them numerically using Octave or R, then we'll come back to the problem.

\subsection{Interlude: calculating the inverse of a matrix using Octave}

The following ad hoc Octave code will create a matrix for $A_B$, take its
inverse, and left multiply it with the full constraint matrix $A$.

\begin{verbatim}
octave:1> AB = [ 1, 1, 0; 2, 1, 1; -1, 1, 0]
AB =

   1   1   0
   2   1   1
  -1   1   0

octave:2> A = [1, 1, 1, 0, 0; 2, 1, 0, 1, 0; -1, 1, 0, 0, 1]
A =

   1   1   1   0   0
   2   1   0   1   0
  -1   1   0   0   1

octave:3> ABinv = inv(AB)
ABinv =

   0.50000   0.00000  -0.50000
   0.50000   0.00000   0.50000
  -1.50000   1.00000   0.50000

octave:4> ABinv*A
ans =

   1.00000   0.00000   0.50000   0.00000  -0.50000
   0.00000   1.00000   0.50000   0.00000   0.50000
   0.00000   0.00000  -1.50000   1.00000   0.50000

\end{verbatim}

\subsection{Interlude: calculating the inverse of a matrix using R}

The following ad hoc R code will create a matrix for $A_B$, take its
inverse, and left multiply it with the full constraint matrix $A$.
Note the use of \verb|%*%| to do the matrix multiplication; simply
using \verb|*| does element-by-element multiplication.  The function
\verb|solve()|, for our purposes here, inverts the matrix.

\begin{verbatim}
> AB = t(matrix(
+ c(1, 1, 0,  2, 1, 1,  -1, 1, 0),
+ nrow = 3,
+ ncol = 3))
> AB
     [,1] [,2] [,3]
[1,]    1    1    0
[2,]    2    1    1
[3,]   -1    1    0
> A = t(matrix(
+ c(1, 1, 1, 0, 0,  2, 1, 0, 1, 0,  -1, 1, 0, 0, 1),
+ nrow=5, ncol=3))
> A
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    1    1    0    0
[2,]    2    1    0    1    0
[3,]   -1    1    0    0    1
> ABinv <- solve(AB)
> ABinv
     [,1] [,2] [,3]
[1,]  0.5    0 -0.5
[2,]  0.5    0  0.5
[3,] -1.5    1  0.5
> ABinv %*% A
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    0  0.5    0 -0.5
[2,]    0    1  0.5    0  0.5
[3,]    0    0 -1.5    1  0.5
\end{verbatim}

\subsection{Now Back to Our Regularly Scheduled Program}

So how did we get that canonical form in
Eqs.~\ref{eq:fav-basis-124-obj} to \ref{eq:fav-basis-124-nonneg}?
There's about 3.5 pages of linear algebra in the textbook...

Fundamentally, take the program (P):
\begin{align}
  & \operatorname{max}\{z(\vec{x}) =  \tvec{c}\vec{x} + z_a\} \\
  \textrm{subject to} & \nonumber\\
  & A\vec{x} = \vec{b} \\
  & \vec{x} \ge \vec{0}
\end{align}
find your next preferred basis $B$, and our new, equivalent program is
\begin{align}
  & \operatorname{max}\{z(\vec{x}) = \tvec{y}\vec{b} + z_a +
  (\tvec{c} - \tvec{y}A)\vec{x}\} \\
  \textrm{subject to} & \nonumber\\
  & A_B^{-1}A\vec{x} = A_B^{-1}\vec{b} \\
  & \vec{x} \ge \vec{0}, \\
  \textrm{where} & \nonumber\\
  & \vec{y} = A_B^{-\intercal}\vec{c_B}.
\end{align}

(Note that $A_B^{-\intercal}$ is the transpose of the inverse of the
matrix $A$, and the inverse operation and transpose operation commute,
so $A_B^{-\intercal} = (A_B^{-1})^\intercal = (A_B^\intercal)^{-1}$.)

Or, I prefer to write it
\begin{align}
  & \operatorname{max}\{z'(\vec{x}) =  \tvec{c'}\vec{x} + z'_a\} \\
  \textrm{subject to} & \nonumber\\
  & A'\vec{x} = \vec{b'} \\
  & \vec{x} \ge \vec{0}
\end{align}
where
\begin{align}
\vec{y} & = A_B^{-\intercal}\vec{c_B} \label{eq:y}\\
\tvec{c'} & = \tvec{c} - \tvec{y}A \\
z'_a & = \tvec{y}\vec{b} + z_a \\
A' & = A_B^{-1}A \label{eq:Aprime}\\
\vec{b'} & = A_B^{-1}\vec{b}.
\end{align}

This, then, is the basic process of canonization (has nothing to do
with the Catholic saints):
\begin{enumerate}
\item rewrite the conditions using the procedure above, to satisfy the
  first condition;
\item rewrite the objective function using $\vec{y}$ from
  Eq.~\ref{eq:y}.
\end{enumerate}

Apply all of this to
Eqs.~\ref{eq:fav-obj-2nd-time}-\ref{eq:fav-nonneg-2nd-time}
and you will get
Eqs.~\ref{eq:fav-basis-124-obj}-\ref{eq:fav-basis-124-nonneg}.


\subsection{Additional Example}

Due to popular request, in class we went through a second example of
creating the canonical form for a different basis.  Still the same problem,

\begin{equation}
  \operatorname{max} z(\vec{x}) = (2,3,0,0,0)(x_1,x_2,x_3,x_4,x_5)^\intercal
  \label{eq:fav-obj-3rd-time}
\end{equation}
subject to
\begin{align}
  \left(
    \begin{array}{ccccc}
      1 & 1 & 1 & 0 & 0 \\
      2 & 1 & 0 & 1 & 0 \\
      -1 & 1 & 0 & 0 & 1 
    \end{array}
  \right)
  \vec{x} =
  \left(
    \begin{array}{c}
      6 \\
      10 \\
      4
    \end{array}
  \right) \\
  x_1, x_2, x_3, x_4, x_5 \ge 0.
  \label{eq:fav-nonneg-3rd-time}
\end{align}

Note that in this case, $z_a = 0$.

As written, this is already in canonical form for $B = \{ 3, 4, 5 \}$.
The requested basis was $B = \{ 2, 4, 5 \}$, which gives

\begin{align}
A_B = 
  \left(
    \begin{array}{ccccc}
      1 & 0 & 0 \\
      1 & 1 & 0 \\
      1 & 0 & 1 
    \end{array}
  \right)
\end{align}

and 

\begin{align}
A_N = 
  \left(
    \begin{array}{ccccc}
      1 & 1 \\
      2 & 0 \\
      -1 & 0 
    \end{array}
  \right)
\end{align}

We are going to need both $A_B^{-1}$ and $A_B^{-\intercal}$, so let's
calculate them:

\begin{verbatim}
octave:1> A = [1, 1, 1, 0, 0; 2, 1, 0, 1, 0; -1, 1, 0, 0, 1]
A =

   1   1   1   0   0
   2   1   0   1   0
  -1   1   0   0   1

octave:2> AB = [1, 0, 0; 1, 1, 0; 1, 0, 1]
AB =

   1   0   0
   1   1   0
   1   0   1

octave:3> ABinv = inv(AB)
ABinv =

   1   0   0
  -1   1   0
  -1  -0   1

octave:4> ABinvT = ABinv'
ABinvT =

   1  -1  -1
   0   1  -0
   0   0   1
\end{verbatim}

(In Octave, the single quote mark, or apostrophe, gives you the
transpose of a matrix or vector.  This is not to be confused with the
``prime'' mark $'$ that I put on variables as I update them to new values.)

In Eq.~\ref{eq:Aprime}, we saw that the rewritten constraint equation uses
$A' = A_B^{-1}A$ and $\vec{b'} = A_B^{-1}\vec{b}$, so let's calculate those:

\begin{verbatim}
octave:5> Aprime = ABinv*A
ans =

   1   1   1   0   0
   1   0  -1   1   0
  -2   0  -1   0   1

octave:6> b = [6;10;4]
b =

    6
   10
    4

octave:7> bprime = ABinv*b
ans =

   6
   4
  -2
\end{verbatim}

That gives us our constraints.  Now, to rewrite the objective function
is a little more complicated.  Recall $\vec{c_B}$ is the subset of
basis elements from the objective function coefficients $\vec{c}$
corresponding to our desired basis $B$, so in this case, elements 2, 4
and 5 of $(2, 3, 0, 0, 0)$, then we calculate $\vec{y} =
A_B^{-\intercal}$,

\begin{verbatim}
octave:8> CB = [3; 0; 0]
CB =

   3
   0
   0


octave:9> y = ABinvT*CB
y =

   3
   0
   0
\end{verbatim}

Now we can calculate the new $c' = c^\intercal - y^{\intercal}A$,

\begin{verbatim}

octave:12> c = [2; 3; 0; 0; 0]
c =

   2
   3
   0
   0
   0


octave:23> cprime = c' - y'*A
cprime =

  -1   0  -3   0   0


\end{verbatim}

Finally, calculate $z'_a = \tvec{y}\vec{b} + z_a$,
\begin{verbatim}
octave:26> y'*b
ans =  18
\end{verbatim}
giving us (hopefully) the revised linear program
\begin{align}
  & \operatorname{max}\{z'(\vec{x}) =  \tvec{c'}\vec{x} + z'_a\} \\
  \textrm{subject to} & \nonumber\\
  & A'\vec{x} = \vec{b'} \\
  & \vec{x} \ge \vec{0}
\end{align}
where
\begin{align}
A' & = \left(
    \begin{array}{ccccc}
   1 &  1 &  1 &  0 &  0 \\
   1 &  0 & -1 &  1 &  0 \\
  -2 &  0 & -1 &  0 &  1
  \end{array}
\right) \\
\vec{b'} & = \left(
    \begin{array}{c}
      6 \\
      4 \\
      -2
    \end{array}
    \right) \\
\vec{c'} & = \left(
    \begin{array}{c}
      -1 \\
      0 \\
      -3 \\
      0 \\
      0
    \end{array}
    \right) \\
z'_a & = 18
\end{align}

If all has gone well, that has changed our problem statement from one
that was canonical in the 3, 4, 5 basis to one that is canonical in
the basis 2, 4, 5.  Confirming our conditions: columns 2, 4, and 5
form an identity matrix ($A'_B = I$) (check!); the corresponding
entries in our objective function are zero ($\vec{c'_B} = \vec{0}$)
(check!).

Whew!  It worked.

\subsection{The Basic Solution}

(This is material not covered in class.)

Every basis is associated with a unique \emph{basic solution}.  The
elements of the basic solution $\vec{x_a}$ related to the basis $B$,
which we will call the vector $\vec{x_B}$, come from the equation
\begin{align}
\vec{x_B} = A_B^{-1}\vec{b}.
\end{align}
The elements not related to $B$ are part of the partition $N$, of
course, and all of those are 0, $\vec{x_N} = \vec{0}$.  Put
$\vec{x_B}$ and $\vec{x_N}$ back together in the proper order to get
$\vec{x_a}$.

\section{The Simplex Algorithm}

\subsection{The Convex Polytope}

A \emph{polytope} is an $n$-dimensional polygon (2-D) or polyhedron
(3-D).  It is \emph{convex} if a line drawn between any two points on
the surface of the polytope passes only through its interior.

A polytope can be created by cutting the space with a set of planes.
In a linear program, each constraint defines a plane, and the total
set defines the polytope.  For our purposes here, the polytope doesn't
have to be bounded in all dimensions and directions.  Before any
constraints, we start with the unbounded polytope that is the positive
quadrant of our space.  If all of the constraints are necessary,
adding the constraints one by one will chop off some of the space and
create a new face for the polytope, keeping it convex and adding some
number of vertices and edges in the process.

Ultimately, the possible solution space of a linear program with a
feasible solution is a polytope.  One (or more) of the vertices of the
polytope will represent the optimal solution.  Thus, we can find that
vertex by starting at any vertex and sliding along the edges until we
find the optimal one.  This sliding is our basic \emph{simplex
  operation}.  (Mathematicians will tell you it's not a precise name,
that ``simplex'' really means the $n$-dimensional triangle, but it's
the common name, so we'll go with it.)


\emph{Question: Does the basic feasible solution we begin with
  \emph{have} to be a vertex?  Why won't any point in the interior of
  the polytope do?  In three dimensions, running in any direction will
  first encounter a face, then a second run along the face will find
  an edge, and the third run along the edge will find a vertex.  Will
  these three steps suffice in 3 dimensions?  Perhaps the basis
  rotation in the core algorithm isn't good enough to make the
  algorithm work right, unless we first find a vector normal to the
  face we encounter on the first run?  What about $n$ dimensions?}


\subsection{Anzen Dai-ichi}

The idea of the simplex algorithm is pretty simple: repeat a simplex
operation until no more simplex operations improve the result, and
you're done.  But there are a couple of things we need to worry about:

\begin{enumerate}
\item The second simplex operation might partly undo the work of the
  first!
\item At each vertex, we have several possible choices for which way
  to go next.  If we choose poorly:
  \begin{enumerate}
  \item we might be inefficient in finishing; or
  \item worse, we might wind up looping forever!
  \end{enumerate}
\item We need to be careful to recognize when we're on an edge that
  extends forever (that is, the solution is unbounded).
\item Finally, of course, we need some way to recognize when we are
  done, and have found the (or an) optimal solution.
\end{enumerate}

\subsection{The Basic Algorithm}

See the separate file uploaded to SFS for pseudocode for the basic
algorithm, slightly reformatted from p. 70 of textbook.

\section{Homework}

See the separate file uploaded to SFS.

\section{References for this Lecture}

\begin{enumerate}
\item
  \verb|https://www.gnu.org/software/octave/doc/interpreter/Three_002dDimensional-Plots.html|
  has examples of the Octave function \verb|patch|.
\item \verb|https://www.gnu.org/software/octave/doc/interpreter/Graphics-Object-Properties.html|
\item \verb|https://jp.mathworks.com/help/matlab/ref/patch.html| is
  help in Japanese on using the \verb|patch| function in Octave.
\item I recommend that you avoid learning more about finding the basic
  solution for right now; it's really messy.  But if you want to, I
  found the video \verb|https://www.youtube.com/watch?v=eYa5aJZhfpg|
  helpful.
\end{enumerate}

\end{document}
