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

\DeclareMathOperator*{\E}{\mathbb{E}}
\let\Pr\relax
\DeclareMathOperator*{\Pr}{\mathbb{P}}

\newcommand{\eps}{\epsilon}
\newcommand{\inprod}[1]{\left\langle #1 \right\rangle}
\newcommand{\R}{\mathbb{R}}

\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}{Scribe: #4}{Lecture #1}}

\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{observation}[theorem]{Observation}
\newtheorem{proposition}[theorem]{Proposition}
\newtheorem{definition}[theorem]{Definition}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{fact}[theorem]{Fact}
\newtheorem{assumption}[theorem]{Assumption}

% 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{22 --- Nov 11, 2014}{Fall 2014}{Prof.\ Eric Price}{Taewan Kim}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                      1. Overview                   %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Overview}

In the last few lectures we covered
\begin{enumerate}
\item Fourier Transform
\item Sparse Fourier Transform
\item Fourier RIP
\end{enumerate}

In this lecture, a new topic 'Oblivious Subspace Embeddings' is covered, especially algorithms introduced by Clarkson and Woodruff \cite{CW13} for regression and low rank approximation problems.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         2. Application
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Application}
Oblivious Subspace Embedding (OSE) is a tool for faster numerical linear algebra. There are two possible applications where OSE can be applied: regression and low rank approximation.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                  2.1 Regression                    %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Regression}
\textbf{Problem Statement} Given $A \in \mathbb{R}^{n\times d}$ and $b\in \mathbb{R}^n$, Find $x\in \mathbb{R}^d$ minimizing $\|Ax-b\|_2$. $(n \gg d)$

$A$ is given data composed of $n$ rows of size $d$ which indicates the $d$ different features. And for those $n$ items, vector $b$ is composed of $n$ outcomes corresponding to each $1 \times d$ feature vector. By finding solution $x\in \mathbb{d}$ minimizing $\|Ax-b\|_2$, we can find an approximate linear mapping between $A$ and $b$ via $x$: $Ax \approx b$.

This problem can be relaxed by allowing $\epsilon$ error:
\begin{equation*}
\text{Find~} x'~s.t.~\|Ax'-b\|_2 \leq (1+\epsilon)\min_x \|Ax-b\|_2
\end{equation*}

An algorithm to find optimal solution of the regression problem ($\epsilon=0$) is using Moore-Penrose pseudoinverse \footnote{Details of Moore-Penrose pseudoinverse can be found in Wikipedia or Chapter 4 of Laub, Alan J. {\em Matrix analysis for scientists and engineers}. Siam, 2005.}.

\textbf{Algorithm}
\begin{enumerate}
\item $x=A^+b$, ($A^+$ is pseudoinverse of $A$) 
\item When $rank(A)=d$ and $d\ll n$, $A^+ = (A^TA)^{-1}A^T$
\end{enumerate}

Time complexity of this algorithm to calculate $x=A^+b$ is $\mathcal{O}(d^2n+d^3)=\mathcal{O}(d^2n)$ when $d\ll n$.
To speed up, sparsity of $A$ can be utilized when $A$ is sparse. If $nnz(A)$ represents the number of nonzero elements in $A$, time complexity can be improved to $\mathcal{O}(d\cdot nnz(A) +d^3)$. The actual computation should be done as follows: compute $A^Tx$ first, which gives $\mathcal{O}(nnz(A))$, then compute $(A^TA)^{-1}$ which gives $\mathcal{O}(d^3)$, and finally compute $(A^TA)^{-1}(A^Tx)$.

However, by using OSE of \cite{CW13} one can achieve:
\begin{itemize}
\item $\mathcal{O}(nnz(A))+\tilde{\mathcal{O}}(d^3/\epsilon^2)$
\item $\mathcal{O}(nnz(A)\log(1/\epsilon))+\tilde{\mathcal{O}}(d^3\log(1/\epsilon))$\\
(Here, $\tilde{\mathcal{O}}(f)\triangleq f\cdot\log^{\mathcal{O}(1)}(f)$)
\end{itemize}


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%                 2.2 Low rank Approximation         %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\subsection{Low Rank Approximation}
\textbf{Problem Statement} Given a matrix $A\in \mathbb{R}^{n\times n}$, find a matrix $B$ with $rank(B)=k$ which minimizes $\|A-B\|_F^2$.

This low rank approximation problem with Frobenius norm can also be relaxed by allowing $\epsilon$ error:
\begin{equation*}
\text{Find~} B'~s.t.~\|A-B'\|_F^2 \leq (1+\epsilon)\min_{\substack{B\\rank(B)=k}} \|A-B\|_F^2
\end{equation*}

When $\epsilon=0$, Singular Value Decomposition (SVD) gives the best rank-$k$ approximation of $A$ by selecting top $k$ singular values and corresponding singular vectors. SVD requires $\mathcal{O}(n^3)$ of computational time.

However by using Power method/subspace iteration:
\begin{itemize}
\item Each iteration takes $\mathcal{O}(n^2k)$ time.
\item For Frobenius norm approximation, bound is not known.
\item By allowing spectral error, $\tilde{\mathcal{O}}(n^2k/\epsilon^2)$ is possible per iteration.
\end{itemize}

Also, utilizing OSE can give better time bound:
\begin{itemize}
\item $\mathcal{O}(nnz(A))+\tilde{\mathcal{O}}(nk^2/\epsilon^4+k^3/\epsilon^5)$
\end{itemize}

For a dense matrix $A$, rank-$k$ matrix approximation using random projection was introduced by \cite{Sarlos06,CW09}.

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%         3. Oblivious Subspace Embedding            %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

\section{Oblivious Subspace Embedding}
\begin{definition}
Defined on parameters $(m,n,d,\epsilon,\delta)$. An Oblivious Subspace Embedding (OSE) is a distribution on matrices $S\in \mathbb{R}^{m \times n}$, s.t. $\forall$ $d$-dimensional subspace $U$ of $\mathbb{R}^n$, with probability $1-\delta$ over $S$, we have $\forall x \in U$ that $\|Sx\|_2=(1\pm \epsilon)\|x\|_2$
\end{definition}

\subsection{Regression with OSE}

Now, we can solve the problem in easier way with lower dimension using OSE. Rather than solving $x^*=\arg\min_x\|Ax-b\|$, solve:
\begin{align*}
x'=&\arg\min_x \|SAx-Sb\|\\
=&\arg\min_x \|S(Ax-b)\|
\end{align*}
where $(Ax-b)\in Col(A\circ b)$. \\
($Col(A\circ b)$ means a column space of $A$ adjoined with the vector $b$)

Then from the definition of OSE,
\begin{equation}
\frac{\|Ax'-b\|}{\|Ax^*-b\|}\leq \left(\frac{1+\epsilon}{1-\epsilon}\right)\frac{\|S(Ax'-b)\|}{\|S(Ax^*-b)\|}\leq \left(\frac{1+\epsilon}{1-\epsilon}\right)\ \lesssim 1+3\epsilon
\end{equation}

Computational time is determined by "$Embedding~time + Solve(m,d)$", where $Solve(m,d)$ represents the time to solve new regression problem with size $m\times d$ of $SA$ and $m\times 1$ vector $Sb$.

One example of OSE is Gaussian random matrix which can be defined as:
\begin{equation}
S_{i,j}=\mathcal{N}(0,1/m)
\end{equation}
With Gaussian OSE, $m=\mathcal{O}(d/\epsilon^2)$. Therefore, embedding requires $\mathcal{O}(mnd)=\mathcal{O}(d^2n/\epsilon^2)$ and $Solve(d,m)$ requires $\mathcal{O}(d^3/\epsilon^2)$ computational time. So, total time is $\mathcal{O}(d^2n/\epsilon^2+d^3/\epsilon^2)$.

\subsection{Fast Johnson-Lindenstrauss}
Now, we introduce an important lemma, which is called Johnson-Lindenstrauss (JL) lemma.
\begin{definition}[Johnson-Lindenstrauss Lemma]
If $m=\mathcal{O}((1/\epsilon^2)\log(1/\delta))$, then

\centering $\forall x$, $\|Ax\|_2^2=(1\pm\epsilon)\|x\|_2^2$\hspace{2em} w.p. $1-\delta$
\end{definition}
Think as this way: given $d$-dim. subspace $U$, take $\epsilon$-net: $C=(1+1/\epsilon)^d$ points. If $m=\mathcal{O}((1/\epsilon^2)\log(1/\delta))$, then all are preserved, i.e. $C$ can be covered.
\begin{align*}
x=x_1+&\epsilon x_2 +\epsilon^2 x_3 + \cdots \text{ for }x_1,\cdots \in C\\
\Rightarrow \|Ax\|_2 &\geq \|Ax_1\|-\epsilon\|Ax_2\|-\epsilon^2\|Ax_3\|-\cdots\\
&\geq 1-\epsilon-(1+\epsilon)(\epsilon+\epsilon^2+\cdots)\\
&\geq 1-3\epsilon
\end{align*}

Faster version of Johnson-Lindenstrauss embedding technique was introduced by \cite{KW11}:

If $A$ has RIP of order $k$, then $AD$ has $(\epsilon,2^{-k})$ JL property, where 
\begin{equation*}
D=\left[ \begin{array}{cccc}
\pm 1& & &\\
& \pm 1 & &\\
& & \ddots &\\
& & & \pm1
\end{array} \right]
\end{equation*}.

Last class, it is shown that $F_{\Omega \in [n]}$ satisfies $(k,\epsilon)$ RIP if $|\Omega|\geq (1/\epsilon^2) k\log^4n$. So, if $m=|\Omega|$ is greater than $(d/\epsilon^2)\log(1/\epsilon) \log^4 n$, then subspace embeddings with $m=(d/\epsilon^2)\log^5n$, and computational time is $n\log n$.
So, with Fast JL, embedding requires $\mathcal{O}(dn\log n)$ and $Solve(m,d)$ requires $\mathcal{O}((d^3/\epsilon^2)\log^5n)$.

\subsection{\cite{CW13}}
To improve the complexity, \cite{CW13} used the sparsity of $A$.
In each column of $S$, exactly one element has $\pm 1$ value defined with hash functions:
\begin{align*}
h&: [n] \rightarrow [m] &\leftarrow 2\text{-independent}\\
\sigma &: [n] \rightarrow \{ \pm 1\} &\leftarrow 4\text{-independent}
\end{align*}
therefore OSE matrix $S$ is defined as,
\begin{equation*}
S_{h(i),i}=\sigma_i
\end{equation*}

Let's prove that above $S$ is OSE by showing:
\begin{equation*}
a,b\in \mathbb{R}^n \Rightarrow \langle Sa,Sb\rangle \approx \langle a,b\rangle
\end{equation*}
 
\begin{proof}
Denote $\delta_{r,i}=\mathbb{I}_{h(i)=r}$ (indicator function).
\begin{align*}
\langle Sa,Sb \rangle &= \sum_{r=1}^{m}\left[\left( \sum_{i=1}^{n}\delta_{r,i}\sigma_{r,i}a_i \right)\left( \sum_{j=1}^{n}\delta_{r,j}\sigma_{r,j}b_i \right)\right]\\
&=\left[\sum_{i=1}^{n}a_i b_i \left( \sum_{r=1}^{m}\delta_{r,i}^2 \sigma_{r,i}^2\right) \right] +\sum_{r=1}^{m}\sum_{i\neq j}\delta_{r,i}\delta_{r,j}\sigma_{r,i}\sigma_{r,j}a_i b_j\\
&=\langle a,b \rangle + \sum_{r=1}^{m}\sum_{i\neq j}\delta_{r,i}\delta_{r,j}\sigma_{r,i}\sigma_{r,j}a_i b_j\\
\Rightarrow \mathbb{E}[\langle Sa,Sb \rangle] &= \langle a,b \rangle
\end{align*}
Now let's consider the variance, $Var[\langle Sa,Sb \rangle ]$. (This proof can be referred to \cite{NN13})
\begin{align*}
(Var[\langle Sa,Sb \rangle ])^2&=\sum_{r=1}^{m}\sum_{i\neq j} \mathbb{E}\left[\sigma_{r,i}^2 \delta_{r,j}^2 (a_i^2 b_j^2 + a_i b_j a_j b_i) \right]\\
&\left( \begin{array}{cl}
\multicolumn{2}{l}{\because \text{Consider } (r,i),(r,j),(r',i'),(r',j')}\\
r=r' \text{ or } \{i,j\}=\{i',j'\}&\rightarrow \mathbb{E}[~\cdot~]\neq 0\\
\text{Otherwise} &\rightarrow \mathbb{E}[~\cdot~]=0 \text{ by independence.}
\end{array} \right)\\
\Rightarrow Var[\langle Sa,Sb \rangle ]&=\frac{1}{m}\sum_{i\neq j} \left( a_i^2 b_j^2 +a_i b_j a_j b_i\right)\\
&\leq \frac{2}{m}\sum_{i\neq j}a_i^2 b_j^2\\
&\leq \frac{2}{m}\sum_{i, j}a_i^2 b_j^2 = \frac{2}{m}\|a\|_2^2 \|b\|_2^2
\end{align*}

Let $U\in \mathbb{R}^{n\times d}$ have orthonormal columns. We want,
\begin{align*}
&\|SUx\|_2 = (1 \pm \epsilon)\|x\|_2 ~~\forall x\in \mathbb{R}^d\\
\Leftrightarrow ~&x^TU^TS^TSUx = (1 \pm \epsilon)x^Tx\\
\Leftrightarrow ~&\|U^TS^TSU-I\|_2 \leq \epsilon\\
\Leftarrow ~&\|U^TS^TSU-I\|_F^2 \leq \epsilon^2
\end{align*}
So it is sufficient to show for Frobenius norm case.
\begin{align*}
(U^TS^TSU)_{i,j} &= \langle SU_i, SU_j \rangle~~(U_i:~i^{th} \text{ column of } U)\\
I_{i,j} &= \langle U_i , U_j \rangle
\end{align*}
Also,
\begin{gather*}
\forall i,j~ \mathbb{E}[(U^TS^TSU-I)_{i,j}^2] \leq \frac{2}{m}\\
\Rightarrow \mathbb{E}[\|U^TS^TSU-I\|_F^2] \leq \frac{2d^2}{m}\leq 2\epsilon^2\\
\Rightarrow \|U^TS^TSU-I\|_2 \leq \epsilon
\end{gather*}
which shows that $\|SUx\|_2 = (1 \pm \epsilon)\|x\|_2 ~~\forall x\in \mathbb{R}^d$, i.e. $S$ is OSE.
\end{proof}
With this setting of $S$ by \cite{CW13}, complexity can be achieved to $\mathcal{O}(
nnz(A)+(d^3/\epsilon^2)\log^5 (d/\epsilon))$, which is $\mathcal{O}(
nnz(A))+\tilde{\mathcal{O}}(d^3/\epsilon^2)$.

Following Table compares the computational time for introduced algorithms when applied to regression problem. ($\mathcal{O}$ notation is omitted.)

\begin{table}[h]
\centering
\begin{tabular}{r c c c}
No OSE: & \multicolumn{3}{c}{ $d\cdot nnz(A)+d^3~/~d^2n+d^3$}\\[.5em]
with OSE: &$Embedding$ &+ &$Solve(d,m)$\\[.5em]
Gaussian: &$mnd=d^2n/\epsilon$ & + &$d^3/\epsilon^2$\\[.5em]
Fast JL: &$dn\log n$ &+ &$(d^3/\epsilon^2)\log^5n$\\[.5em]
C-W: &$nnz(A)$ &+ &$Solve(d,d^2/\epsilon^2)$\\
& & & $=d^4/\epsilon^2 \leftarrow$ bad!\\
& & & $\rightarrow(d^3/\epsilon^2)\log^5 (d/\epsilon)$
\end{tabular}
\caption{Comparing complexities for various algorithms for regression}
\end{table}

\newpage
\bibliographystyle{alpha}

\begin{thebibliography}{42}

\bibitem[CW09]{CW09}
Clarkson, Kenneth L., and David P. Woodruff.
\newblock Numerical linear algebra in the streaming model.
\newblock {\em Proceedings of the forty-first annual ACM symposium on Theory of computing}, ACM, 2009.

\bibitem[CW13]{CW13}
Clarkson, Kenneth L., and David P. Woodruff.
\newblock Low rank approximation and regression in input sparsity time.
\newblock {\em Proceedings of the forty-fifth annual ACM symposium on Theory of computing}, ACM, 2013.

\bibitem[KW11]{KW11}
Krahmer, Felix, and Rachel Ward.
\newblock New and improved Johnson-Lindenstrauss embeddings via the restricted isometry property.
\newblock {\em SIAM Journal on Mathematical Analysis} 43.3 (2011): 1269-1281.

\bibitem[NN13]{NN13}
Nelson, Jelani, and Huy L. Nguyên.
\newblock OSNAP: Faster numerical linear algebra algorithms via sparser subspace embeddings.
\newblock {\em Foundations of Computer Science (FOCS), 2013 IEEE 54th Annual Symposium on}. IEEE, 2013.

\bibitem[Sarlos06]{Sarlos06}
Sarlos, Tamas. 
\newblock Improved approximation algorithms for large matrices via random projections.
\newblock {\em Foundations of Computer Science, 2006. FOCS'06. 47th Annual IEEE Symposium on}, IEEE, 2006.

\end{thebibliography}

\end{document}