\documentclass[11pt]{article}
\usepackage{amsmath,amssymb,amsthm,mathtools}

\DeclareMathOperator*{\E}{\mathrm{E}}
\DeclareMathOperator*{\Var}{\mathrm{Var}}
\DeclareMathOperator*{\argmin}{arg\,min}
\DeclareMathOperator*{\median}{\mathrm{median}}
\DeclareMathOperator*{\mean}{\mathrm{mean}}

\newcommand{\eps}{\epsilon}
\newcommand{\inprod}[1]{\left\langle #1 \right\rangle}
\newcommand{\bbR}{\mathbb{R}}
\newcommand{\R}{\mathbb{R}}
\newcommand{\poly}{{\rm poly}}
\newcommand{\norm}[2]{||#1||_{#2}}

\newcommand{\handout}[5]{
  \noindent
  \begin{center}
  \framebox{
    \vbox{
      \hbox to 5.78in { {\bf CS 395T: Sublinear Algorithms } \hfill #2 }
      \vspace{4mm}
      \hbox to 5.78in { {\Large \hfill #5  \hfill} }
      \vspace{2mm}
      \hbox to 5.78in { {\em #3 \hfill #4} }
    }
  }
  \end{center}
  \vspace*{4mm}
}

\newcommand{\lecture}[4]{\handout{#1}{#2}{#3}{Scribes: #4}{Lecture #1}}


\newcommand{\brk}{\vspace{ 2 mm }}
\newcommand{\Brk}{\vspace{ 4 mm }}

\newtheorem{thm}{Theorem}[section]
\newtheorem{claim}[thm]{Claim}
\newtheorem{lemma}[thm]{Lemma}
\newtheorem{define}[thm]{Definition}
\newtheorem{cor}[thm]{Corollary}
\newtheorem{obs}[thm]{Observation}
\newtheorem{remark}[thm]{Remark}
\newtheorem{example}[thm]{Example}
\newtheorem{construction}[thm]{Construction}
\newtheorem{conjecture}{Conjecture}
\newtheorem{THM}{Theorem}
\newtheorem{COR}[THM]{Corollary}
\newtheorem{question}{Question}
\newtheorem{fact}[thm]{Fact}
\newtheorem{prop}[thm]{Proposition}
% 1-inch margins, from fullpage.sty by H.Partl, Version 2, Dec. 15, 1988.
\topmargin 0pt
\advance \topmargin by -\headheight
\advance \topmargin by -\headsep
\textheight 8.9in
\oddsidemargin 0pt
\evensidemargin \oddsidemargin
\marginparwidth 0.5in
\textwidth 6.5in

\parindent 0in
\parskip 1.5ex

\begin{document}

\lecture{14 --- October 13, 2016}{Fall 2016}{Prof.\ Eric Price}{Cody Freitag, William Swartworth }

\section{Compressed Sensing}

In compressed sensing, the overall goal is to observe an object using as few measurements as possible.  In today's lecture, we'll motivate compressed sensing with an application to observing images. Then we'll give the precise setup for compressed sensing.  Lastly we'll give an interative algorithm to solve compressed sensing.

\subsection{Single-Pixel Cameras}

We will consider taking a picture as a motivating example.  A normal camera uses $1000 \times 1000$ pixels with 1 photosensitive element per pixel. This is typically cheap for normal images since we can use silicon.  However, silicon doesn't work if you want to capture infrared images, making it really expensive.  As a result, we may try to decrease the number of photosensitive elements needed at the expense of taking a greater number of measurements per element.  

An image can be represented as a vector $x = (x_1, \ldots, x_n)$ where entries of $x$ correspond to individual pixels.  Instead of getting a measurement for each pixel, we could instead get $\sum_i b_i x_i = \langle b ,  x \rangle$ for some $b \in \{0,1\}^n$ (where we think of the $b_i$'s as lens filters on the camera).  We could do this many times with different filters and about $\approx n/2$ ones per filter $b$.  Then with around $n$ different $b$, we can reconstruct $x$ just by inverting the $b$ matrix.  So we would only need one photosensitive element, but we'd need $n$ measurements now.

This leads to the fundamental problem in compressed sensing: choose a matrix $A$ such that observing $Ax$ allows $x$ to be recovered.

\subsection{Questions to Consider}

Do we really need $n$ measurements?

The answer is yes if we require exact recovery of every possible $x$.  Many situations are more relaxed than this.  For example it might be possible to recover $x$ approximately from $Ax$, even if $Ax$ is much smaller than $x$, if $x$ is known to have some structure.  In fact, if $x$ is an image, then we expect to be able to compress it anyway (using JPEG for example).  So it's reasonable to hope that $Ax$ could be substantially smaller than $x$.

\Brk

How do we model compressibility?

We can often think of compressible data as being sparse with respect to some basis.  This is a familiar idea in the context of lossy compression.  To compress an image we start by changing basis---a typical choice is the wavelet basis.  If we've chosen a good basis, then many entries will be negligibly small, and we approximate our image as a sparse vector in our chosen basis.


\subsection{Compressed Sensing - Specifying the Problem}

We begin by choosing $A \in \bbR^{m \times n}$ where generally $m$ is much smaller than $n$.  

Making our measurements gives us some vector $y = Ax + e$, where we assume that $\norm{x}{0} \le k$ (where the $0$-norm is defined to be the sparsity).  The vector $e$ should be thought of as an error term which models imprecise measuring equipment.

Our goal is to recover $x$ from $y$. In other words, we would like to find $\widehat x \approx x$ given $y$ (and knowledge of the $A$ that we chose).

This can be summarized as follows:
\begin{align*}
& \text{Choose } A \in \bbR^{m \times n} \\ 
& \text{Observe } y = Ax + e, \text{ where } \norm{x}{0} \le k \text{ and } e \text{ is noise} \\
& \text{Compute } \widehat x \approx x \text{ from } (A,y)
\end{align*}

\section{Algorithms}

Many algorithms exist for compressed sensing:
\begin{itemize}
\item Convex Optimization

\begin{itemize}
\item L1 minimization, LASSO, Dantzig selector
\end{itemize}

\item Iterative methods

\begin{itemize}
\item CoSAMP, AMP, OMP, IHT
\end{itemize}
\end{itemize}

Many of these algorithms use a matrix $A$ satisfying RIP($O(k)$, $0.1$) to find $\widehat x$ with error $\norm{\widehat x - x}{2} \lesssim \norm{e}{2}$. Fortunately we can construct RIP($O(k)$, $0.1$) in many ways:
\begin{itemize}
\item random gaussian with $m = O(k \log (n/k) )$

\item any JL matrix for $\delta = 2^{-O(k \log (n/k) )}$

\item matrices with low coherence [columns are $a_1, \ldots, a_n$ and $\frac{| \langle a_i, a_j \rangle |}{\sqrt{\norm{a_i}{} \cdot \norm{a_j}{}}} < \frac{1}{k}$ for all $i, j$.  This is easy to check but requires $m > k^2$.]

\item random rows of Fourier, $O(k \log n \log^2 k)$

\item explicity $k^{2-\epsilon}$ row matrices exist.

\end{itemize}

We'll look at IHT (Iterative Hard Thresholding) today.

\subsection{IHT}

First suppose that $A$ satisfies the $(3k,\epsilon)$-RIP.

This means that for all $3k$-sparse $x$ in $\R^n$ we have 
\begin{align*}
\norm{Ax}{}^2 &= (1+ \epsilon) \norm{x}{2}^2\\
x^T A^T A x &= (1 + \epsilon) x^T x\\
|x^T (A^T A - I) x| &\le \epsilon x^T x.\\
\end{align*}
In terms of the operator norm, 
$$\norm{(A^T A - I)_{S \times S}}{} \le \epsilon, \forall S \subset [n], |S| \le 3k.$$

So the RIP implies that $A^T A$ is approximately the identity on small sets of coordinates.

\Brk

In the compressed sensing problem, we're given $Y = A x + e$, where $x$ is $3k$-sparse, and $\norm{e}{}$ is ``small''.

Let $z = A^T y = A^T Ax + A^T e$.  Then $A^T e$ should be small noise, and $A^T A x$ should be approximately $x$ as long as $x$ is sparse.  So we expect that $||z-x||$ is small.

To be more precise, let $S$ be any set of size at most $3k$, which we assume contains the support of $x$. Then we have 
\begin{align*}
\norm{(z-x)_S}{} 
	&\le \norm{(A^TA - I)_{S \times S} \cdot x_S}{} + \norm{(A^Te)_S}{} \\
	&\le \epsilon \cdot \norm{x}{} + \norm{A^T_{S \times [m]}}{} \cdot \norm{e}{}.
\end{align*}

The bound on $\norm{(A^Te)_S}{}$ was shown in the next class.

\Brk

\begin{lemma}  Let $x,z \in \bbR^n$, with $x$ $k$-sparse and with support $H$.  Let $S$ be the set of indices corresponding to the top $k$ elements of $z$.  Then $\norm{x - z_S}{2}^2 \le 5 \cdot \norm{(x-z)_{H \cup S}}{2}^2$.
\end{lemma}

\begin{proof} Pair up $i \in H \setminus S$ and $j \in S \setminus H$.  We know $|z_j| > |z_i|$, so we just need $x_i^2 + z_j^2 \le 5 \cdot ((x - z_i^2) + z_j^2)$.  There are two cases.  Either $|z_i| > |x_i| / 2$, which implies $x_i^2 + z_j^2 \le 4z_i^2 + z_j^2 \le 5 z_j^2$.  Or, $|z_i| \le |x_i| / 2$, which implies $x_i^2 + z_j^2 \le 4 (x_i - z_i)^2 + z_j^2$.
\end{proof}

Now for $z = A^T y = A^T (Ax) + A^T e$, let $T$ denote the top $k$ elements of $z$. Applying the lemma to $T \cup {\rm supp}(x)$ we get
\begin{align*}
\norm{x - z_T}{2} 
	&\le \sqrt{5} \cdot \norm{(x-z)_{T \cup S}}{2}  \\
	&\le \sqrt{5} ( \epsilon \norm{x}{} + 3 \norm{e}{}).
\end{align*} 
For $\epsilon = 0.1$, this is at most $\norm{x}{} / 4 + 3 \norm{e}{}$.

Now let $y^\prime = A(x - z_T) + e = y^\prime - Az_T$.

The same analysis as before says $z^\prime = A^T y^\prime = A^T A (x - z_T) + a^T e$ has 
\begin{align*}
\norm{(z^\prime - (x - z_T))+S}{} 
	&= \norm{(x - (z_T + z^\prime))_S}{} \\
	&\le \epsilon \norm{x- z_T}{} + (1 + \epsilon) \norm{e}{}.
\end{align*}

Set $x^{(2)}$ to be the top $k$ elements of $z_T + z^\prime$.  According to the lemma, 
\begin{align*}
\norm{x - x^{(2)}}{} 
	&\le \sqrt{5} \norm{x - (z_T + z^\prime)_{{\rm supp}(x^{(2)}) \cup {\rm supp}(x)} }{} \\
	&\le \sqrt{5} \epsilon \norm{x - z_T}{} + (1+\epsilon)\norm{e}{} \\
	&\le 1/4 \cdot (\norm{x - x^{(1)}}{} + 3 \norm{e}{}).
\end{align*}

Now we just repeat to a get a sequence of successively better recoveries of x.  So set $$z^{(t+1)} = A^T(y - Ax^{(t)}),$$ and $x^{(t+1)} =$ top $k$ entries of $z^{(t+1)}$.  Our analysis gives us the following:
$$\norm{x^{(t+1)} - x} \le \frac{1}{4} \norm{x^{(t)} - x}{} + 3 \norm{e}{}.$$

Therefore if $\norm{x^{(t+1)} - x}{} > 12 \norm{e}{}$, then $\norm{x^{(t+1)} - x}{} \le \frac{1}{2} \norm{x^{(t)} - x}{}.$  Hence, after $I = \log (\norm{x}{} / \norm{e}{})$ iterations, we get a recovery $\widehat x = x^{(I)}$ with $\norm{\widehat x - x}{} \le 12 \norm{e}{}$.  

The running time of our algorithm is $$O\left(\log \frac{\norm{x}{} }{ \norm{e}{}} \cdot \text{matrix/vector multiply}\right).$$


\end{document}