
\documentclass{article}
\usepackage[margin=1in]{geometry}
\usepackage{subfig}
\usepackage{amssymb}
\usepackage{amsmath,amsthm}
\usepackage{graphicx,hyperref}

\newcounter{lecnum}
\newcommand{\R}{\mathbb{R}}
\newcommand{\transpose}{\mathtt{T}}   

\newcommand{\lecture}[4]{
   \newpage
   \setcounter{lecnum}{#1}
   \noindent

   \begin{center}
   \framebox{
      \vbox{\vspace{2mm}
    \hbox to 16cm { {\bf CS602 Applied Algorithms
                        \hfill 2019-20 Sem II} }
       \vspace{4mm}
       \hbox to 16cm { {\Large \hfill Lecture #1 -  II: #2  \hfill} }
       \vspace{2mm}
       \hbox to 16cm { {\it Scribe: #4  \hfill  Lecturer: #3} }
      \vspace{2mm}}
   }
   \end{center}
   \vspace*{4mm}
}


\newtheorem{theorem}{Theorem}[lecnum]
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{definition}[theorem]{Definition}




\begin{document}

%\lecture{**LECTURE-NUMBER**}{**DATE**}{**LECTURER**}{**SCRIBE**}
\lecture{18}{March 24}{Rohit Gurjar}{Priyansh Kimtee}

Now, we will see the ellipsoid method and an interior point method to solve linear programs.

\section{Ellipsoid Method}
 The main advantage of this method over other two methods is that it does not require an explicit description of the Linear Program(LP), 
 we just need a separation oracle.
 
 \subsection{Separation Oracle}
Let  $C \subseteq \mathbb{R}^n$ be a closed convex set and $x \in \mathbb{R}^n$ not in $C$. 
Then $x$ and $C$ can be strictly separated by a \textbf{separating hyperplane}.(As seen in earlier lectures).

Separation oracle either says that the given point is inside the convex body or it if outside then gives a certificate for this fact, which is a separating hyperplane.

\subsection{Ellipsoid}
A 2D ellipse centered at origin can be described as:
\[\frac{x_1^2}{b_1^2}+\frac{x_2^2}{b_2^2}\leq{1}\]
where  \(b_1\) and \(b_2\) are the lengths of major and minor axes.

An ellipsoid is a higher dimensional generalization of an ellipse. 
An $n$ dimensional ellipsoid centered at origin can be described as :
\[\frac{x_1^2}{b_1^2}+\frac{x_2^2}{b_2^2}+.....+\frac{x_n^2}{b_n^2}\leq{1}\]
where  $b_1, b_2, \dots, b_n$ are the lengths of the axes.
To get an ellipsoid in a general orientation, one can apply a rotation transformation on the coordinates. 

\subsection{Feasibility}
As we already know that Feasibility and Optimization are equivalent, we will just describe how the ellipsoid method is used to find a feasible point
for the LP, i.e., a point inside the polytope.
% whether the given LP is feasible or not i.e the polytope is empty or non-empty.
\paragraph{Algorithm}
\textbf{Invariant:} An ellipsoid $E$  that  contains the given polytope.  Let $z$ denote the the center of the current ellipsoid. 
\begin{enumerate}
    \item Initially, we construct an ellipsoid such that the polytope is 
    completely inside it. This can be done by taking 
    a sphere centered at the origin with large enough radius. 
  Initially $z= (0,0, \dots, 0)$.
    \item We give the point $z$ and the polytope (convex body) to the separation oracle. 
    
    \textbf{if} the point is inside the polytope \textbf{then} done\\
    \textbf{else} it gives a separating hyperplane, say $(a,b)$ such that $a^\top z > b$, but $a^\top x < b $ for each $x$ in the polytope.
    
    \item Next we find a new ellipsoid which contains that half of the current ellipsoid where the polytope lies i.e., $  \{x \in E \mid a^\top x < b \}$. 
%    To elaborate, if the polytope lies $< $ side of the hyperplane (i.e., every point in the polytope satisfies $a^\top x < b$)
%That is, we want to find a new ellipsoid that contains the half-ellipsoid 
   % also the new ellipsoid passes through the point of intersection of the separating hyperplane and the current 
%    ellipsoid. Now the new point P is the center of the new ellipsoid.
    
    \item Go to step 2 with the new ellipsoid $E$ and its center $z$.
\end{enumerate}

\subsubsection{Running time}
\begin{itemize}
    \item If the initial conditions are good like the polytope should be in the initial sphere (which can be done with the help of constraints and estimating the distance of the farthest point from the origin in the polytope). %and taking that as the radius of the sphere).
    \item Volume of the polytope should not be too small (reason shown in section \ref{analysis})
    \item If the above condition are satisfied then the algorithm runs in polynomial time.
\end{itemize}

\subsubsection{Analysis of ellipsoid algorithm}\label{analysis}
We need to show that the volume of the ellipsoid shrinks by a constant factor in each round. 
That is,
\[\frac{Vol(E_{i+1})}{Vol(E_i)}\leq{\delta} \]
for some $\delta <1$.
For any ellipsoid \(E_{i}\) and polytope P we can say  \(E_i\supseteq P\) therefore \[Vol(P)\leq{Vol(E_i)}\].
If the algorithm takes $r$ rounds to finish then 
\[\delta^r \geq \frac{Vol (P)}{Vol (E_0)}.\]
taking log on both the sides we get
\[r \leq \log_\frac{1}{\delta} \left( \frac{Vol (E_0)}{Vol (P)} \right).\]
%where r denotes the number of times we need to use the Separation oracle.

Thus, the runtime is small when the volume of the polytope is large enough.
In particular,  we should have $Vol(P) \neq 0$, i.e., the polytope should be full dimensional. 

\textbf{Note:} If the LP is explicit, the separation oracle is very easy, we just have to check each constraint for a point one by one and if the point is not feasible then one of the constraint give you the separating hyperplane. 
For many combinatorial problems, it turns out that there is an exponentially large number of constraints in the linear program.
Ellipsoid method can still be used to optimize if a separation oracle can be designed. 
\\
\\
\textbf{Points to think:}
\begin{itemize}
    \item How to construct new ellipsoid which satisfy the required conditions.
    \item How to bound the volume of the ellipsoid in each iteration and also the volume of polytope.
\end{itemize}


\section{Interior point method}
The starting point of the interior point methods is the Newton Raphson Method.
It also uses the duality of convex programs. 
\subsection{Newton-Raphson Method}
Newton-Raphson method is used to find the roots of a polynomial function.

Consider a polynomial function $P(x)$,
\begin{itemize}
    \item Start with a arbitrary point say \(x_0\) and draw a tangent at this point on the curve $P(x)$. This tangent can be thought of as 
    a linear approximation of $P(x)$ at $x_0$ %which gives us the approximate value at \(x_0\).
    \item Next let the point where tangent intersect x-axis be \(x_1\) (i.e., root of the linear approximation).
        If $P^\prime(x)$ is the derivative of $P(x)$ then we can write % the slope be denoted by \(P^\prime(x_0)\) at point \(x_0\)
        \[P^\prime(x_0) = \frac{P(x_0)}{x_0-x_1}\]
        \[x_1-x_0=-\frac{P(x_0)}{P^\prime(x_0)}\]
        for any general \({(i+1)}^{th}\) step,
        \[x_{i+1}-x_{i}=-\frac{P(x_i)}{P^\prime(x_i)}\]
     \end{itemize}
   
        \begin{figure}
            \centering
            \includegraphics{download.png}
            \caption{Newton-Raphson Method  (Image source: \url{ http://fourier.eng.hmc.edu/e176/lectures/NM/node20.html} ) }
            \label{fig:nrm}
        \end{figure}

It is possible that this algorithm may not converge if starting point is not chosen correctly;
it might get stuck in a cycle.
So, the starting point needs to be chosen carefully. 
We are not going into the details of how to choose the staring point. 

\paragraph{Multivariate Newton-Raphson:}
The Newton-Raphson method can be generalized to multivariate and multi-valued polynomials. 
Let  
$P:\mathbb{R}^n \rightarrow \mathbb{R}^m$ be such a polynomial and we are looking for its root. 
In other words, we want to solve the following system of polynomial equations. 

\begin{equation}
    \begin{aligned}
    P_1(x_1,x_2, \dots, x_n)&=0\\
    P_2(x_1,x_2, \dots, x_n)&=0\\
    &\vdots\\
    P_m(x_1,x_2, \dots, x_n)&=0
    \end{aligned}
\end{equation}

To run Newton-Raphson method, we need the multi-variate analogue of the derivative,
which is given by the  following Jacobian matrix $J(x)$. 


\[
\begin{bmatrix}
  \frac{\partial P_1}{\partial x_1} & 
    \frac{\partial P_1}{\partial x_2} & 
     \cdots &
    \frac{\partial P_1}{\partial x_n} \\[1ex] % <-- 1ex more space between rows of matrix
  \frac{\partial P_2}{\partial x_1} & 
    \frac{\partial P_2}{\partial x_2} & 
     \cdots&
    \frac{\partial P_2}{\partial x_n} \\[1ex]
   && \vdots &\\
  \frac{\partial P_m}{\partial x_1} & 
    \frac{\partial P_m}{\partial x_2} & 
     \cdots&
    \frac{\partial P_m}{\partial x_n}
\end{bmatrix}
\]

%Therefore we use the jacobian matrix to get the \(P^\prime(x)\),
Similar to the univariate Newton-Raphson, we can write the following relation between 
the $i^{th}$ and $(i+1)^{th}$ point. Let $J= J(x=x_i)$.
\begin{eqnarray}
& J  (x_{i+1}-x_i)=-P(x_i) \nonumber \\ 
\implies & x_{i+1}-x_{i}=-J^{-1} P(x_i) \label{eq:newton-step}
\end{eqnarray}
If the matrix \(J\) is not invertible then we take its pseudo inverse \(J^+\)
where 
\begin{equation*}
\begin{aligned}
& J^+=(J^TJ)^{-1} J^T\\
\therefore\hspace{1cm}
& x_{i+1}-x_i=-J^+ P(x_i)
\end{aligned}
\end{equation*}
%
Now, suppose we want to use Newton-Raphson to solve the following LP. 
%We cannot immediately apply the above as there are also inequalities. 
\begin{eqnarray}
\min \quad c^\intercal x \nonumber \\
Ax=b \label{eq:primal1}\\
x \geq  0 \label{eq:primal2}
\end{eqnarray}
Its dual LP would be:
\begin{eqnarray}
\max \quad b^\intercal y \nonumber \\
A^\intercal y \leq c \label{eq:dual}
\end{eqnarray}
Recall that the complementary slackness is given as
%
\begin{equation}
x_i (a_{i} ^\intercal y-c_i)=0 \text{ for each } 1 \leq i \leq n \label{eq:comp},
\end{equation}
%
where $a_i$ is the $i$-th column of $A$.
First step is to do away with the \emph{minimization} question. 
We know that any pair $(x,y)$ of feasible solutions that also satisfy (\ref{eq:comp}) 
will be optimal solutions. 
So, we can ignore the minimization/maximization question and just ask for a solution 
that simultaneously satisfies (\ref{eq:primal1}), (\ref{eq:primal2}), (\ref{eq:dual}), and (\ref{eq:comp}).
Note that (\ref{eq:comp}) is not linear but a degree-2 polynomial equation.
And so Newton-Raphson comes into the picture.
However, we cannot immediately apply Newton-Raphson for this as (\ref{eq:primal2}) and (\ref{eq:dual})
are inequalities.
%If we are able to solve equation (5),(6),(8),(9)
%then we can find (x,y),but equation (6) is 
The idea is to eliminate the inequality (\ref{eq:primal2}) 
% by introducing a \textbf{Barrier function} .
and  introduce a \textbf{Barrier Function} that helps us put the inequality (\ref{eq:primal2})  implicitly  in the objective function. 
%
We can do this by creating a function that greatly increases the objective if a constraint 
is about to be violated. 
%

Our new objective function will be $f(x) = c^\intercal x - \mu \sum_{i}\ln{x_i}$. 
%
Note that  $-\ln{x_i}$ is increases when $x_i$ decreases, and in fact,
the rate of increase goes higher and higher as $x_i$ goes close to zero. 
%
So, if we are in the quadrant $x > 0$ and are  trying to minimize $f(x)$, we are forced to be away from 
the hyperplanes $x_i=0$. 
%
However, by introducing this barrier function, we have changed our objective function. 
%
The parameter $\mu$ gives us a tradeoff: when $\mu$ is large, more weight is given to satisfying
the inequality (\ref{eq:primal1}); when $\mu$ is small then $f(x)$ is close to being the actual objective function.
%
% so due to barrier function we can drop the constraint \(x \geq 0\) also when \(\mu=0\), then
%a minimum to the barrier function will be a minimum to our original LP.
Our new program becomes
\begin{eqnarray}
\min \quad f(x) = c^\intercal x - \mu \sum_{i}\ln{x_i}  \nonumber \\ 
\text{ subject to } 
 Ax = b. \label{eq:primal}
\end{eqnarray}
We can write the dual of this program similar to how we write it for linear programs. 
%
For LP, the $c$ vector in the objective goes in the dual constraints. Here, instead of $c$,
the partial derivative vector of the objective function will go in the dual constraints. Note that when $\mu=0$,
the derivative is just the vector $c$.
The dual program is
\begin{eqnarray*}
\max \quad  b^\intercal y \\
\text{ subject to }
 A^\intercal y=\nabla f
\end{eqnarray*}

where \(\nabla {f}\) is the partial derivative vector of $f$ w.r.t.\ \(x_1,x_2, \dots, x_n\).
That is, 
\[
\nabla f = ( \dots, c_i-\frac{\mu}{x_i}, \dots)
\]
So, the dual program can be written as 
\begin{eqnarray}
&\max \quad  b^\intercal y& \nonumber \\
\text{ subject to }&
 a_i^\intercal y= c_i-\frac{\mu}{x_i} &\text{ for } 1\leq i \leq n \label{eq:dual2}.
\end{eqnarray}

Again to get rid of the minimization/maximization question, we can use complementary slackness conditions,
which is actually empty (as there are no non-negativity constraints on $x$ and $y$).
%
So, if we can just find a pair  $(x,y)$ that satisfies (\ref{eq:primal}) and (\ref{eq:dual2})
simultaneously, they will be optimal solutions for the above programs. 
Hence, we need to solve the following system of polynomial equations.     
\begin{eqnarray}
     & \mu=x_i(c_i-a_{i} ^ \intercal y) & \text{ for } 1\leq i \leq n \label{eq:1}\\
      &Ax=b   \label{eq:2} &
\end{eqnarray}
%Now this is a polynomial equation and now 
Now, we can apply Newton-Raphson method to find a solution for the above system.
%\textbf{Note:}Here we do not have complementary slackness condition as we do not have non-negativity constraints on x and y variables.
%We start with some $x>0$ and a large value  and $y$ satisfying
 We start with large value of \(\mu\) as first we want to give more weightage to satisfying $x>0$.
 %
 With each Newton-Raphson step, we will slowly decrease $\mu$, as small $\mu$ means our objective function
 is close to the actual objective. 
The following summarizes the algorithm: Fix $\rho < 1$.
\begin{enumerate}
\item Start with some   $x,y$ (as discussed before, starting point needs to be chosen carefully).
\item Take a Newton-Raphson step (take Jacobian of (\ref{eq:1}) and (\ref{eq:2}) and apply rule (\ref{eq:newton-step})). 
\item Update $\mu \leftarrow \rho*\mu$ and go to 2. 
\end{enumerate} 

\end{document}