\documentclass[final,ps,azure]{prosper}
\usepackage{algorithm}
%\usepackage[section]{algorithm}
\usepackage{algorithmic}
\usepackage{amsmath}
\usepackage{amssymb} 
\usepackage{pstcol} % To use the standard `color' package with Seminar
\usepackage{semcolor}
\usepackage{eufrak}
\usepackage{epsfig}
%\usepackage{epsfig}
\def\halmos{\hfill $\Box$\\ }
\newcommand{\abs}{{\rm abs\,}}
\newcommand{\vc}{{\rm vec\,}}
\newcommand{\Diag}{{\rm Diag\,}}
\newcommand{\diag}{{\rm diag\,}}
\newcommand{\tr}{{\rm trace }}
\newcommand{\trace}{{\rm trace\,}}
\newcommand{\Zero}[1]{{\bf 0}_{#1 \times #1}}
\newcommand{\Eye}[1]{{\bf I}_{#1 \times #1}}
\def\RE{\mathbb R}
\def\Rn {\RE^n}
\def\Rnp {\RE^n_+}
\def\Rm {\RE^m}
\def\Rnn {\RE^{n\times n}}
\def\Rmn {\RE^{m\times n}}
\def\X {\mathbb X}
\def\SE {\mathbb S}
\def\ME {\mathbb M}
\def\Sn {\SE^n}
\def\Mn {\ME^n}
\def\Snp {\SE^n_+}
\def\A {{\cal A}}
\newcommand{\sri}{\mbox{\rm sri\,}}  
\newcommand{\gra}{\mbox{\rm gra\,}}  
\newcommand{\cone}{\mbox{\rm cone\,}}  
\newcommand{\closu}{\mbox{\rm cl\,}}  
\newcommand{\clconv}{\overline{\mathrm{conv}}\,}
\newcommand{\range}{\mbox{\rm range\,}}  
\newcommand{\rank}{\mbox{\rm rank\,}}  
\newcommand{\fixx}{\mbox{\rm Fix\,}}  
\newcommand{\D}[2]{[{\EuFrak D}#1(#2)]}
\newcommand{\bigD}[2]{\bigl[{\EuFrak D}#1(#2)\bigr]}
\newcommand{\Df}{\D{f}{x}}
\newcommand{\inner}[2]{{\langle #1 , #2 \rangle}}

\slideCaption{{\textcolor{green}{University of Waterloo Tutte Seminar - C\&O}}}

\title{\textcolor{blue}{Projections Methods \\and \\Linear Programs}}
\subtitle{(Reflection-Projections)}
\author{\Large Serge Kruk \\ (Joint with Heinz Bauschke of Guelph)}
\email{\large sgkruk@acm.org}
\institution{\Large Oakland University}
\newcommand{\partials}[2]{\frac{\partial #1}{\partial #2}}


\begin{document}
\maketitle
\begin{slide}{Overview}
  \begin{itemize}
  \item Feasibility Problem
  \item History
  \item New Algorithm(s)
  \item Theoretical convergence
  \item Implementation
  \item Experimentation
  \item Open Questions
  \end{itemize}
\end{slide}

\begin{slide}{Feasibility Problem}
  Given convex sets:
  \begin{displaymath}
    C_1, C_2, \ldots, C_n 
  \end{displaymath}
  
  \begin{itemize}
  \item (Easy?) Find 
    \begin{displaymath}
      x \in C_1 \cap C_2 \cap \ldots \cap C_n
    \end{displaymath}
  \item (Harder?) Or certify no such $x$ exists
  \end{itemize}
\end{slide}

%% \begin{slide}{Projection Algorithms}
%%   Given \\
%%   \quad sets
%%   \begin{displaymath}
%%     C_1, C_2, \ldots, C_n \subset \X
%%   \end{displaymath}
%%   \quad projectors
%%   \begin{displaymath}
%%     P_1, P_2, \ldots, P_n
%%   \end{displaymath}
%%   \quad starting point
%%   \begin{displaymath}
%%     x^0 \in \X
%%   \end{displaymath}
%%   Iterate
%%   \begin{displaymath}
%%     \boxed{
%%       x^{k+1} := x^k + \alpha \sum\omega_i(P_i(x^k)- x^k)
%%     }
%%   \end{displaymath}
%% \end{slide}

\begin{slide}{Projections\\ (von Neuman 1933)}
  Given affine spaces
  \begin{displaymath}
    C_1, \ldots, C_n \quad\mbox{where}\quad 
    C:= C_1\cap\ldots\cap C_n \ne \emptyset
  \end{displaymath}
  \quad any
  \begin{displaymath}
    x^0 \in \X
  \end{displaymath}
  The sequence
  \begin{displaymath}
    x^{k+1} := P_{(k\bmod n) + 1}(x^k)
  \end{displaymath}
  converges to
  \begin{displaymath}
    P_C(x^0)
  \end{displaymath}
\end{slide}

\begin{slide}{Alternating Projections}
    \epsfig{file=altprog.eps,height=3in} 
\end{slide}

\begin{slide}{Reflections \\(Motzkin-Schoenberg 1954)}
  To solve
  \begin{displaymath}
    \inner{a_i}{x}\le b_i \qquad 1\le i \le m
  \end{displaymath}
  use partial reflections
  \begin{displaymath}
        x^{k+1} := \lambda R_{(k\bmod m) + 1}(x^k) \qquad 0<\lambda <2
  \end{displaymath}
  to converge to a solution
\end{slide}

\begin{slide}{Applications}
  \begin{itemize}
  \item Image reconstruction
  \item Tomography
  \item Filter design
  \item \ldots
  \item Dual cone LP
  \end{itemize}
  Common: large problem / low accuracy
\end{slide}

\begin{slide}{Classical application}
  \begin{displaymath}
    \min\{\inner{c}{x} \mid Ax=b, x\ge 0\}
  \end{displaymath}

  Phase I simplex method
  \begin{displaymath}
    x\in \{x\mid Ax=b, x\ge 0\}
  \end{displaymath}
  (and $x$ basic)
\end{slide}

\begin{slide}{IP Application}
  Phase one for feasible interior-point algorithms?
  \begin{displaymath}
    \min \inner{x}{x}
  \end{displaymath}
  Subject to
  \begin{eqnarray*}
    \inner{A_i}{x} &=& b_i, \qquad i\in[1, \ldots m]\\
    x &\in& K
  \end{eqnarray*}
  where $K$ is a closed, convex, pointed, cone.
\end{slide}


\begin{slide}{Why that application?}
  \begin{itemize}
  \item Small problems
    \begin{itemize}
    \item Feasible algorithm
    \item No need for self-dual embedding
    \end{itemize}
  \item Large problem
    \begin{itemize}
    \item Primal-only (or dual-only)
    \item Work in the nullspace of $A$
    \end{itemize}
  \end{itemize}
\end{slide}

\begin{slide}{The final goal}
  \begin{displaymath}
    \boxed{\min\{\inner{c}{x} \mid \inner{A_i}{x} = b_i, x\in K\}}
  \end{displaymath}
  \vspace{1cm}
  \begin{itemize}
  \item Find a feasible point
  \item Express sparse basis for nullspace of $A_i$
  \item Collapse problem to subspace
  \item Solve by feasible interior-point
  \end{itemize}
\end{slide}

\begin{slide}{The final goal}
  \begin{displaymath}
    \boxed{\min\{\inner{c}{x} \mid \inner{A_i}{x} = b_i, x\in K\}}
  \end{displaymath}
  \vspace{1cm}
  \begin{itemize}
  \item \textcolor{red}{Find a feasible point}
  \item Express sparse basis for nullspace of $A_i$
  \item Collapse problem to subspace
  \item Solve by feasible interior-point
  \end{itemize}
\end{slide}


\begin{slide}{Additional structure}
  Let
  \begin{displaymath}
    C := K \cap C_1 \cap \cdots \cap C_N
  \end{displaymath}
  where $K$ is
  \begin{itemize}
  \item Convex
  \item Closed
  \item Pointed
  \item Obtuse
  \item Cone
  \end{itemize}
\end{slide}


\begin{slide}{Cone $K$ Obtuse (Goffin)}
\emph{{Positive Polar Cone}}:
\begin{equation*}
  K^\oplus := \{x \in \X : \inf \langle x,K \rangle \geq 0\}
\end{equation*}
\emph{{Negative Polar Cone}}:
\begin{equation*}
  K^\ominus := \{x \in \X : \inf \langle x,K \rangle \le 0\}
\end{equation*}
\vspace{1cm}

\begin{displaymath}
  \boxed{K \mbox{obtuse} \Leftrightarrow   K^\oplus  \subseteq K}
\end{displaymath}
\end{slide}

\begin{slide}{Obtuse cones}
  \emph{{self-dual cones} } are obtuse
  \vspace{1cm}
  \begin{itemize} 
  \item nonnegative orthant $\Rnp \subset \Rn$;
  \item positive semidefinite matrices $\Snp \subset \Sn$
  \item Lorentz cone
  \end{itemize}
  \vspace{.3in}
  Also the ice-cream cones, \ldots
  \begin{displaymath}
    \mbox{ice}(\alpha) := \{(x,r) \in \X \times \RE \mid \|x\| \le \alpha r\}
  \end{displaymath}
\end{slide}

\begin{slide}{What does obtuseness buy?}
    \epsfig{file=reflection.eps,height=3in} 
\end{slide}

\begin{slide}{Overview}
  \begin{itemize}
  \item Feasibility problem
  \item History
  \end{itemize}
\end{slide}

\begin{slide}{Reflection-Projection ($C\ne \emptyset$)}
  \vspace{1cm}
  \begin{algorithmic}
    \STATE $k:=0$\hfill\COMMENT{Iteration index} 
    \STATE Given $x_k \in
    \X$\hfill\COMMENT{Starting point} 
    \WHILE{$\max \big\{d(x_k,K),d(x_k,C_1),\ldots,d(x_k,C_N)\big\}>0$} 
    \STATE $x_k^{++} := R_{K}x_k$\hfill\COMMENT{Reflect into the cone} 
    \STATE $x_{k+1} := P_N P_{N-1}\ldots P_1x_k^{++}$ 
    \hfill\COMMENT{Project cyclically} 
    \STATE $k:=k+1$\hfill\COMMENT{Next iterate} 
    \ENDWHILE
  \end{algorithmic}
\end{slide}
\begin{slide}{Reflection-Projection ($C\ne \emptyset$)}
  \vspace{1cm}
  \begin{algorithmic}
    \STATE $k:=0$\hfill\COMMENT{Iteration index} 
    \STATE Given $x_k \in
    \X$\hfill\COMMENT{Starting point} 
    \WHILE{$\max \big\{d(x_k,K),d(x_k,C_1),\ldots,d(x_k,C_N)\big\}>0$} 
    \STATE $x_k^{++} := \textcolor{red}{R_{K}x_k}$\hfill\COMMENT{Reflect into the cone} 
    \STATE $x_{k+1} := P_N P_{N-1}\ldots P_1x_k^{++}$ 
    \hfill\COMMENT{Project cyclically} 
    \STATE $k:=k+1$\hfill\COMMENT{Next iterate} 
    \ENDWHILE
  \end{algorithmic}
\end{slide}

\begin{slide}{Reflection-Projection ($C\ne \emptyset$)}
  \vspace{1cm}
  \begin{algorithmic}
    \STATE $k:=0$\hfill\COMMENT{Iteration index} 
    \STATE Given $x_k \in
    \X$\hfill\COMMENT{Starting point} 
    \WHILE{$\max \big\{d(x_k,K),d(x_k,C_1),\ldots,d(x_k,C_N)\big\}>0$} 
    \STATE $x_k^{++} := R_{K}x_k$\hfill\COMMENT{Reflect into the cone} 
    \STATE $x_{k+1} := \textcolor{red}{P_N P_{N-1}\ldots P_1x_k^{++}}$ 
    \hfill\COMMENT{Project cyclically} 
    \STATE $k:=k+1$\hfill\COMMENT{Next iterate} 
    \ENDWHILE
  \end{algorithmic}
\end{slide}

\begin{slide}{Convergence Result}

The sequence
\begin{displaymath}
\boxed{
  (y_i):=   \big(x_0,x_0^{++},x_1,x_1^{++},\ldots, x_k, \ldots\big)
}
\end{displaymath}

\textcolor{blue}{converges to some point in $C$}

Heinz Bauschke and SK, 2002
\end{slide}

\begin{slide}{Projector $P_C$}
  Given convex set $C$ in space $\X$
  
  \begin{displaymath}
    P_C : \X \rightarrow C \qquad \|x-P_C(x)\| = d(x,C)
  \end{displaymath}

  Characterisation
  \begin{displaymath}
    P_C(x) \in C \qquad \inner{S-P_C(x)}{x-P_C(x)}\le 0
  \end{displaymath}

  Firmly nonexpansive
  \begin{displaymath}
    \|P_Cx - P_Cy\|^2 + \|(I-P_C)x - (I-P_C)y\|^2 \leq \|x-y\|^2
  \end{displaymath}
  
\end{slide}

\begin{slide}{Reflector $R_K$}
\begin{itemize} 
\item $R_K$ is \emph{{nonexpansive}}: 
\begin{equation*}
\|R_Kx-R_Ky\| \leq \|x-y\|.
\end{equation*}
\item $K$ is obtuse $\Leftrightarrow$ $R(\X) = K$.
\end{itemize}

Note: $R_K(x) = (2P_K - I)(x)$
\end{slide}


\begin{slide}{Fej\'er monotonicity}
Sequence $y_n$
\begin{equation*}
(\forall n \geq 0)(\forall s \in C) \quad \|y_{n+1}-s\| \leq \|y_n-s\|.
\end{equation*}
Then: \hfill
\begin{enumerate} 
\item $(y_n)$ is bounded, and $\big(d(y_n,C)\big)$ converges.
\item $(P_Cy_n)$ converges to some point $\bar{s} \in C$.
\item $(y_n)$ converges to $\bar{s}$ $\Leftrightarrow$\\
all cluster points of $(y_n)$ belong to $C$. 
\end{enumerate}
\end{slide}

\begin{slide}{Proof of Convergence}
  \begin{eqnarray*}
  (y_i)&:=&   \big(x_0,x_0^{++},x_1,x_1^{++},\ldots, x_k, \ldots\big)\\
    C &:=& K \cap C_1 \cap \cdots \cap C_N
  \end{eqnarray*}


  \begin{itemize}
  \item $(y_k)$ is Fej\'er monotone w.r.t. $C$
  \item $(\forall i)$ $d(x_{k}^{++},C_i) \rightarrow 0$
  \item Each cluster point of $(x_{k}^{++})$ lies in $C$
  \item $(x_k^{++})$ converges to some point $\bar{c} \in C$
  \item $(y_k)$ converges to $\bar{c} \in C$
  \end{itemize}
\end{slide}

\begin{slide}{$(y_k)$ is Fej\'er monotone w.r.t. $C$}

On the one hand, $(y_k)$ is generated by cyclic application of
\emph{nonexpansive} operators $T$ drawn from 
$\{R_K,P_1,\ldots,P_N\}$.

On the other hand, $\fixx R_K = K$ and each $\fixx P_i = C_i$.

Altogether, if $c \in C = \fixx R_K \cap \bigcap_{i} \fixx P_i$ and 
$y_{k+1} = Ty_k$, then
\begin{equation*}
  \|y_{k+1}-c\| = \|Ty_k-Tc\| \leq \|y_k -c\|.
\end{equation*}
\vfill

\end{slide}

\begin{slide}{$(\forall i)$ $d(x_{k}^{++},C_i) \rightarrow 0$}

$K$ obtuse $\Rightarrow$  $(x_{k}^{++}) \in K$.

$P_1$ firmly nonexpansive $\Rightarrow$
\begin{equation*}
  \begin{split}
    &\qquad d^2(x_k^{++},C)\\
    &=\|x_{k}^{++}-P_Cx_k^{++}\|^2 \\
    &\geq \|P_1x_k^{++}-P_Cx_k^{++}\|^2 + \|x_k^{++}-P_1x_k^{++}\|^2\\
    &\geq d^2(P_1x_k^{++},C) + d^2(x_k^{++},C_1).
  \end{split}
\end{equation*}
This, \texttt{Step~1} and the Fej\'er Lemma yield
\begin{equation*}
d^2(x_k^{++},C_1) \leq d^2(x_k^{++},C) - d^2(P_1x_k^{++},C) \rightarrow 0.
\end{equation*}
Argue similarly (and inductively) for $i \geq 2$.


\end{slide}


\begin{slide}{Cluster points of $(x_{k}^{++})$ lie in $C$}

In view of \texttt{Step 2} and the closedness of $K$,
the entire sequence $(x_k^{++})$ lies in $K$ and so
do its cluster points

On the other hand, fix $i$. 
Each $d(x_{k}^{++},C_i) \rightarrow 0$, by 
\texttt{Step~2}. Since $d(\cdot,C_i)$ is continuous, each cluster
point of $(x_{k}^{++})$ must belong to $C_i$.

Altogether, each cluster point of $(x_{k}^{++})$ belongs to
$K \cap C_1 \cap \cdots \cap C_N = C$. 

\end{slide}

\begin{slide}{$(x_k^{++}) \rightarrow \bar{c} \in C$}

On the one hand, by \texttt{Step~1}, $(x_k^{++})$ is Fej\'er monotone
w.r.t. $C$. 

On the other hand, by \texttt{Step~3}, all cluster points
of $(x_k^{++})$ belong to $C$.

Altogether, by the Fej\'er Lemma, $(x_k^{++})$ converges to some $\bar{c}$
in $C$. 

\end{slide}

\begin{slide}{$(y_k)$ converges to $\bar{c} \in C$}

  By \texttt{Step~4}, $(x_{k}^{++}) = (y_{k(N+1)+1})$ converges to 
  $\bar{c} \in C$. Hence, since $P_1$ is continuous and 
  $\bar{c} \in C \subset C_1 = \fixx P_1$, we have
  \begin{equation*}
    P_1x_k^{++} = y_{k(N+1)+2} \rightarrow P_1\bar{c}=\bar{c}.
  \end{equation*}
  Continuing in this fashion yields the result. $\Box$
  
\end{slide}

\begin{slide}{Overview}
  \begin{itemize}
  \item Feasibility problem
  \item History
  \item Abstract algorithm
  \end{itemize}
  
\end{slide}

\begin{slide}{Implementation}
  \vspace{1cm}
  Find
  \begin{displaymath}
    x\in K \cap L
  \end{displaymath}

  where
  \begin{displaymath}
   K = \Rnp \qquad \mbox{or} \qquad K = \Snp
  \end{displaymath}
  and
  \begin{displaymath}
    L = \{ x \in \X : \inner{A_i}{x}=b_i, i=1,\ldots, m \}
  \end{displaymath} 
\end{slide}

\begin{slide}{Implementing Projector $P_K$}

$\boxed{K=\Rnp}$ 

\begin{displaymath}
x=[x]_{j=1}^n \Rightarrow [P_K(x)]_j = x_j^+ = \max \{x_j,0\} 
\end{displaymath}
\end{slide}

\begin{slide}{Implementing Projector $P_K$}

$\boxed{K=\Rnp}$ 

\begin{displaymath}
x=[x]_{j=1}^n \Rightarrow [P_K(x)]_j = x_j^+ = \max \{x_j,0\} 
\end{displaymath}

$\boxed{K=\Snp}$ 

\begin{displaymath}
X = \sum_{j} \lambda_{j}u_j u_j^* \Rightarrow 
P_K(X) = \sum_{j:\lambda_j > 0} \lambda_j u_j u_j^*
\end{displaymath}
\end{slide}

\begin{slide}{Implementing Projector $P_K$}

$\boxed{K=\Rnp}$ 

\begin{displaymath}
x=[x]_{j=1}^n \Rightarrow [P_K(x)]_j = x_j^+ = \max \{x_j,0\} 
\end{displaymath}

$\boxed{K=\Snp}$ 

\begin{displaymath}
X = \sum_{j} \lambda_{j}u_j u_j^* \Rightarrow 
P_K(X) = \sum_{j:\lambda_j > 0} \lambda_j u_j u_j^*
\end{displaymath}

\textcolor{red}{Lanczos iterations}

\end{slide}

\begin{slide}{Implementing Reflector $R_K$}

$\boxed{K=\Rnp}$ 

\begin{displaymath}
x=[x]_{j=1}^n \Rightarrow (R_Kx)_j = |x_j|
\end{displaymath}

$\boxed{K=\Snp}$

\begin{displaymath}
X = \sum_{j} \lambda_{j}u_j u_j^* \Rightarrow 
R_K(X) = \sum_{j} |\lambda_j|  u_j u_j^*.
\end{displaymath}

\end{slide}

\begin{slide}{Implementing Projector $P_L$}

\begin{eqnarray*}
  L &=& \{x \in \X \mid\inner{A_i}{x}=b_i, i=1,\ldots,m\}\\
  &=&\{x\in\X\mid \A x=b\},
\end{eqnarray*}

\begin{equation*}
P_L(x) = x - \A^\dagger(\A(x)-b),
\end{equation*}

$\A^\dagger$ is the \emph{{Moore-Penrose inverse}}  of $\A$.
\end{slide}

\begin{slide}{Implementing $P_L$ (Theory)}

$\boxed{\X=\Rn}$ 

\begin{displaymath}
  A^* = QR \Rightarrow A^\dagger =  QR^{-*}
\end{displaymath}
Therefore
\begin{displaymath}
  P_L(x) = x - A^\dagger(Ax-b) = x - QR^{-*}(Ax-b)  
\end{displaymath}
\end{slide}

\begin{slide}{Implementing $P_L$ (Practice)}
  \vspace{2cm}
  \begin{itemize}
  \item Permute rows of $\A$ for stability 
  \item Rank detection
  \item Remove redundant constraints
  \item Detect infeasibility ($L = \emptyset$)
  \end{itemize}
\end{slide}

\begin{slide}{Implementing $P_L$ (Practice)}
  
\begin{equation*}
  A^* =
  \begin{bmatrix}
    Q & Q_0
  \end{bmatrix}
  \begin{bmatrix}
    R & D \\ 0 & 0
  \end{bmatrix} P^*
\end{equation*}

\begin{equation*}
  \begin{split}
    A^\dagger &= A^*(AA^*)^\dagger \\
    &= Q
    \begin{bmatrix}
      R & D
    \end{bmatrix}
    P^* \Bigg[ P
    \begin{bmatrix}
      R^* \\ D^*
    \end{bmatrix}
    Q^*Q
    \begin{bmatrix}
      R & D
    \end{bmatrix}
    P^*
    \Bigg]^\dagger\\
    &= Q
    \begin{bmatrix}
      R & D
    \end{bmatrix}
    P^* \Bigg[ P
    \begin{bmatrix}
      R^*R & R^*D \\ D^*R & D^*D
    \end{bmatrix}
    P^* \Bigg]^\dagger.
  \end{split}
\end{equation*}

\end{slide}

\begin{slide}{Algo $\{x\in\Rn \mid Ax=b\} \cap \Rnp$}
  \begin{algorithmic}
    \STATE Factor $
    \begin{bmatrix}
      A & b
    \end{bmatrix}^* P = QR$
    \STATE $I_{Ab} := (|\diag(R)| > 10^{-13}\|A\|);$
    \STATE Factor $A^*P = QR$
    \STATE $I := (|\diag(R)| >
    10^{-13}\|A\|);$
    \IF{$|I_{Ab}|\ne |I|$} 
    \STATE Quit
    \ELSE
    \STATE $A := P(I)^*A;\; b:=P(I)^*b;$
    \STATE $Q := Q(I); \;R := R(I);\;$
    \ENDIF 
  \end{algorithmic}
  \hfill\Large\textcolor{blue}{\ldots}
\end{slide}
\begin{slide}{\ldots $\{x\in\Rn \mid Ax=b\} \cap \Rnp$}
  \begin{algorithmic}
    \STATE $r := Ax-b$;\; 
    \STATE Solve $R^* y=r;\;x^+:=x-Qy;\;d_l:=\|x-x^+\|;\;$
    \STATE $v := 2x+e;\; a:=b:=0;\;$
    \WHILE{$(\|d_l\|>\mbox{tol} \wedge\frac{\|v-(a-b)\|}{(1+\|v\|)}>\mbox{tol})$} 
    \STATE $v := a-b;\;b := x ;\; x := x^+ ;\; a := x;\;$ 
    \STATE $x := \abs(x^+);$
    \STATE $r := Ax-b$;\; 
    \STATE Solve $R^* y = r;\;x^+ := x - Qy;\;d_l := \|x-x^+\|;$
    \ENDWHILE
  \end{algorithmic}
\end{slide}

\begin{slide}{\ldots $\{x\in\Rn \mid Ax=b\} \cap \Rnp$}
  \begin{algorithmic}
    \STATE $r := Ax-b$;\; 
    \STATE Solve $R^* y=r;\;x^+:=x-Qy;\;d_l:=\|x-x^+\|;\;$
    \STATE $v := 2x+e;\; a:=b:=0;\;$
    \WHILE{$(\|d_l\|>\mbox{tol}\wedge\textcolor{red}{\frac{\|v-(b-a)\|}{(1+\|v\|)}>\mbox{tol}})$} 
    \STATE $v := b-a;\;b := x ;\; x := x^+ ;\; a := x;\;$ 
    \STATE $x := \abs(x^+);$
    \STATE $r := Ax-b$;\; 
    \STATE Solve $R^* y = r;\;x^+ := x - Qy;\;d_l := \|x-x^+\|;$
    \ENDWHILE
  \end{algorithmic}
\end{slide}


\begin{slide}{Implementing $P_L$}

$\boxed{\X=\Sn}$ 

\begin{displaymath}
  \vc : \Sn \rightarrow \RE^{t(n)}
\end{displaymath}

\begin{displaymath}
t(n) = 1 + 2 + \cdots + n = \frac{n(n+1)}{2}
\end{displaymath}

\begin{eqnarray*}
  \inner{A_i}{X} &=& \trace A_iX \\
  &=& [\vc(A_i)]^t \vc(X) \\
  &=& \inner{\vc(A_i)}{\vc(X)}
\end{eqnarray*}
\end{slide}

\begin{slide}{Implementing $P_L$}
  \vspace{2cm}
  \begin{itemize}
  \item Translate ``$\A(x)=b$'' from $\Sn$ to $\RE^{t(n)}$
  \item Handle the translated constraints as before
  \end{itemize}
\end{slide}

\begin{slide}{Overview}
  \begin{itemize}
  \item Feasibility problem
  \item History
  \item Abstract algorithm
  \item Experiments (and dreams)
  \end{itemize}
\end{slide}

\begin{slide}{Experimental setup}
  \begin{itemize}
  \item $\X = \RE^{64}$;
  \item 1,000 random problems 
    \begin{equation*}
      Ax=b, \quad x \geq 0,
    \end{equation*}
    all feasible; 
  \item $2 \leq m \leq 64$
  \end{itemize}
\end{slide}

\begin{slide}{Experiment}

\begin{equation*}
x \mapsto \big( (1-\alpha_L)I + \alpha_L P_L\big) \circ
\big( (1-\alpha_K)I + \alpha_K P_K\big),
\end{equation*}
where
\begin{equation*}
0<\alpha_K \leq 2 \quad \text{and} 
\quad 0<\alpha_L < 2.
\end{equation*}

\begin{itemize} 
\item $\alpha_K=1$ and $\alpha_L=1$ $\Rightarrow$ cyclic projections;
\item $\alpha_K=2$ and $\alpha_L=1$ $\Rightarrow$ reflection-projection.
\end{itemize}

\end{slide}

\begin{slide}{Results (Optimal step)}
\begin{center}
\epsfig{file=contour64.eps,width=4in}
\end{center}

\end{slide}

\begin{slide}{Results (Optimal step)}
  \begin{center}
    \epsfig{file=surfc64.eps,width=4in}
  \end{center}
\end{slide}


\begin{slide}{Iteration count $\alpha_L=1$}
  \epsfig{file=middle.ps,angle=-90,width=4in}
\end{slide}

\begin{slide}{Infeasible start (Traditional)}
  \vspace{.5in}
  \fbox{Simplex phase I} $\rightarrow$ \fbox{Simplex phase II} 
\end{slide}

\begin{slide}{Infeasible start (Interior-point)}
  \vspace{.5in}
  \fbox{Infeasible IP algo} or \fbox{Self-dual-embedding}
\end{slide}

\begin{slide}{Phase I for IP?}
  \vspace{.5in}
  \fbox{Reflection-Projection} $\rightarrow$ 
  \fbox{Feasible IP algo}
\end{slide}

\begin{slide}{Phase I for IP?}
  \vspace{.5in}
  \fbox{Reflection-Projection} $\rightarrow$ 
  \fbox{Feasible IP algo}\\
  \vspace{.5in}
  \begin{itemize}
  \item In theory: \textcolor{blue}{bad idea}
  \item In practice ?
  \end{itemize}
\end{slide}

\begin{slide}{Non-polynomial}
    \epsfig{file=badreflection.eps,height=3in} 
\end{slide}

\begin{slide}{Epelman-Freund (2002)}
  \vspace{1in}
  Some projection methods \textcolor{blue}{are} polynomial
\end{slide}

\begin{slide}{Infeasible start vs Feasible start}
  \vspace{.5in}
  \begin{itemize}
  \item Generate random LP
  \item Solve via infeasible start algo
  \item Find a starting point via Reflection-Projection
  \item Re-solve via feasible start algo
  \end{itemize}
\end{slide}

\begin{slide}{Starting point for primal-dual}
  \vspace{.5in}
  \begin{align*}
    \min&\{\inner{c}{x} \mid \A(x)= b, x\in K\}\\
    \max&\{\inner{b}{y} \mid \A^*(y)+z-C = 0, z\in K^*\}
  \end{align*}
  \end{slide}

  \begin{slide}{Perturbed Optimality}
    \begin{align*}
      \A(x) = b\\
      \A^*(y) +z = c\\
      \inner{c}{x} - \inner{b}{y} = \mu\\
      x \in K\\
      z \in K^*
    \end{align*}
    
  \end{slide}
  \begin{slide}{Central-path conditions}
  \vspace{.5in}
  \begin{displaymath}
    \begin{bmatrix}
      \A & 0 & 0 \\
      c^* & -b ^* & 0 \\
      0 & \A^* & I 
    \end{bmatrix}
    \begin{bmatrix}
      x \\ y \\ z
    \end{bmatrix}
    =
    \begin{bmatrix}
      b \\ \mu \\ c
    \end{bmatrix}
  \end{displaymath}
  \begin{displaymath}
    (x,y,z) \in (K \times \RE^m \times K^*)
  \end{displaymath}
  \begin{itemize}
  \item Easy extension of our proof
  \item Numerical algebra cost is small
  \end{itemize}
\end{slide}

\begin{slide}{Infeasible start vs Feasible start}
  \begin{table}[htbp]
    \centering
    \begin{tabular}{ccc}
      Infeasible  & R-P & Feasible\\ \hline
       11  & 249 &     10\\
       11  & 264 &     10\\
       10  & 249 &       9\\
       10  & 302 &       8\\
       11  & 188 &       9\\
       10  & 234 &       9\\
       10  & 332 &       8\\
       10  & 319 &       8\\
       10  & 248 &       9
    \end{tabular}
    \caption{Iteration counts}
  \end{table}
\end{slide}

\begin{slide}{Solve the lp?}
Why not?
  \begin{displaymath}
    \begin{bmatrix}
      \A & 0 & 0 \\
      c^* & -b ^* & 0 \\
      0 & \A^* & I 
    \end{bmatrix}
    \begin{bmatrix}
      x \\ y \\ z
    \end{bmatrix}
    =
    \begin{bmatrix}
      b \\ c \\ 0
    \end{bmatrix}
  \end{displaymath}
  \begin{displaymath}
    (x,y,z) \in (K \times \RE^m \times K^*)
  \end{displaymath}
\end{slide}

\begin{slide}{Why not?}
  \vspace{.5in}
  Because it takes \textcolor{red}{forever...}
  
  \textcolor{blue}{Condition number of the system?}
 \end{slide}

\begin{slide}{Why not?}
  \vspace{.5in}
  Because it takes \textcolor{red}{forever...}
  
  \textcolor{blue}{Condition number of the system?}
  \begin{itemize}
  \item $F := \{(\A,b) \mid \A(x)=b, x\in K\}$ is feasible
  \item $\rho := \inf\{\|\Delta\A,\Delta b\| \mid
    (\A+\Delta\A,b+\Delta b) \in \partial F\}$
  \item $\kappa := \frac{\|\A,b\|}{\rho}$
  \end{itemize}
\end{slide}


\begin{slide}{Overview}
  \begin{itemize}
  \item Feasibility problem
  \item History
  \item Abstract algorithm
  \item Experiments (and dreams)
  \item Loose ends
  \end{itemize}
\end{slide}


\begin{slide}{Feasibility Problem}
  Given convex sets:
  \begin{displaymath}
    C_1, C_2, \ldots, C_n 
  \end{displaymath}
  
  \begin{itemize}
  \item (Easy?) Find 
    \begin{displaymath}
      x \in C_1 \cap C_2 \cap \ldots \cap C_n
    \end{displaymath}
  \item (Harder?) Or certify no such $x$ exists
  \end{itemize}
\end{slide}

 \begin{slide}{Inconsistent case}
   \vspace{.5in}
 \begin{itemize} 
 \item $A$ and $B$ are two closed convex sets in $\X$; 
 \item projectors $P_A$ and $P_B$; 
 \item reflectors $R_A=2P_A-I$ and $R_B=2P_B-I$. 
 \end{itemize}

 \end{slide}

 \begin{slide}{Gap Vector}

 \begin{equation*}
 v := P_{\closu(B-A)}(0),
 \end{equation*}
 and set
 \begin{equation*}
 E :=  A \cap (B-v) \text{~and~} F :=  B \cap (A+v).
 \end{equation*}
 Then $E = \fixx (P_AP_B)$ and $F = \fixx (P_BP_A)$.\\
 If $e \in E$, then $P_B(e) = P_F(e) = e+v$. \\
 (Similarly for $f \in F$.)
 \end{slide}

 \begin{slide}{The geometry of $A \cap B = \emptyset$}
   \epsfig{file=figure_4.eps,height=3in}
 \end{slide}
 
 \begin{slide}{R-P on infeasible problems}
   
   Iterates are
   \begin{equation*}
     b_{k} = R_Ba_k, \quad a_{k+1} := P_Ab_k,
   \end{equation*}
   
 Then:

 $\bullet$ $\fixx(P_AR_B)=E$ and $\fixx(R_BP_A)=F+v$;\\
 $\bullet$ $(\forall a \in A)(\forall k \geq 0)$:
 \begin{equation*}
 \begin{split}
 \|a_k-a\|^2 \geq& \|a_{k+1}-P_AR_Ba\|^2 + \\
 &\|(b_k-a_{k+1})-(R_Ba-P_AR_Ba)\|^2.
 \end{split}
 \end{equation*}

 \end{slide}


 \begin{slide}{If gap realized}
 \begin{itemize} 
 \item[(i)] $(a_k)$ is Fej\'er monotone w.r.t. $E$; in fact,
 \begin{equation*}
 \|a_k-e\|^2 \geq \|a_{k+1}-e\|^2 + \|(b_k-a_{k+1})-2v\|^2,
 \end{equation*}
 for every $e \in E$ and $k \geq 0$. 
 \item[(ii)] $b_k-a_{k+1} \rightarrow 2v$.
 \item[(iii)] $\big((a_k+a_{k+1})/2\big)$ clusters in $E$.
 \item[(iv)] $(P_Ba_k)$ clusters in $F$. 
 \end{itemize}
 \end{slide}

 \begin{slide}{If gap realized}
 \begin{itemize} 
 \item $(a_k)$ is Fej\'er monotone w.r.t. $E$; in fact,
 \begin{equation*}
 \|a_k-e\|^2 \geq \|a_{k+1}-e\|^2 + \|(b_k-a_{k+1})-2v\|^2,
 \end{equation*}
 for every $e \in E$ and $k \geq 0$
 \item \textcolor{red}{$b_k-a_{k+1} \rightarrow 2v$}
 \item $\big((a_k+a_{k+1})/2\big)$ clusters in $E$
 \item $(P_Ba_k)$ clusters in $F$ 
 \end{itemize}

 $\Rnp$ works but $\Snp$ may not!

 \vfill
 \end{slide}

\begin{slide}{\ldots $\{x\in\Rn \mid Ax=b\} \cap \Rnp$}
  \begin{algorithmic}
    \WHILE{$(\|d_l\|>\mbox{tol} \wedge \textcolor{red}{\frac{\|v-(b-a)\|}{(1+\|v\|)}>\mbox{tol}})$} 
    \STATE $v := b-a;\;b := x ;\; x := x^+ ;\; a := x;\;$ 
    \STATE $x := \abs(x^+);$
    \STATE $r := Ax-b$;\; 
    \STATE Solve $R^* y = r;\;x^+ := x - Qy;\;d_l := \|x-x^+\|;$
    \ENDWHILE
  \end{algorithmic}
  \vspace{.3in}
  Cerfificate of infeasibility
\end{slide}


\begin{slide}{Conclusions}
\vspace{2cm}
\begin{itemize} 
\item Family of projection methods
\item Reflection-Projection - Obtuse Cone
\item Better than cyclic projections
\item Inconsistent case 
\item IP solver  ``phase one''
\item Relation to condition number
\item Adaptive step length ($> 2$)
\end{itemize}

\end{slide}

\begin{slide}{References}

$\bullet$ H. H. Bauschke and S. G. Kruk.\\
\emph{The Method of Reflection-Projection for 
Convex Feasibility Problems with an Obtuse Cone},
March 2002, submitted. 
Available at 
\begin{center}
\verb+www.optimization-online.org+
\end{center}
$\bullet$ \texttt{GNU Octave} and \texttt{Matlab} files available at
\begin{center}
\verb+www.oakland.edu/~kruk/research+
\end{center}

\end{slide}


\end{document}

%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "tutte"
%%% End: