\documentclass[11pt]{article}
\usepackage{latexsym}
\usepackage{amsmath,amssymb,amsthm}
\usepackage{epsfig}
\usepackage{enumerate}
\usepackage[right=0.8in, top=1in, bottom=1.2in, left=0.8in]{geometry}
\usepackage{setspace}
\spacing{1.06}

\newcommand{\handout}[5]{
  \noindent
  \begin{center}
  \framebox{
    \vbox{\vspace{0.25cm}
      \hbox to 5.78in { {COMS E6998-9:\hspace{0.12cm}Algorithmic
          Techniques for Massive Data} \hfill #2 }
      \vspace{0.48cm}
      \hbox to 5.78in { {\Large \hfill #5  \hfill} }
      \vspace{0.42cm}
      \hbox to 5.78in { {#3 \hfill #4} }\vspace{0.25cm}
    }
  }
  \end{center}
  \vspace*{4mm}
}
\newcommand{\lecture}[4]{\handout{#1}{#2}{#3}{Scribes:\hspace{0.08cm}#4}{Lecture #1}}

\theoremstyle{plain}
\newtheorem{theorem}{Theorem}
\newtheorem{corollary}[theorem]{Corollary}
\newtheorem{lemma}[theorem]{Lemma}
\newtheorem{observation}[theorem]{Observation}
\newtheorem{example}[theorem]{Example}
\newtheorem{claim}[theorem]{Claim}
\newtheorem{fact}[theorem]{Fact}
\newtheorem{assumption}[theorem]{Assumption}

\theoremstyle{definition}
\newtheorem{definition}{Definition}

\theoremstyle{remark}
\newtheorem{recall}{Recall}
\newtheorem{remark}{Remark}
\newtheorem{note}{Note}

\newcommand{\E}{\mathbb{E}}
\newcommand{\var}{\text{Var}}

\DeclareMathOperator*{\argmin}{argmin}

\begin{document}

\lecture{11 -- Applications of Dimension Reduction}{Oct 22, 2015}{Instructor:\hspace{0.08cm}\emph{Alex Andoni}}{\emph{Marshall Ball}}

Today we looked at two applications of dimension reduction for improving the time complexity of two classical problems  in {\sc P}:
\begin{enumerate}[(a)]
\item Matrix Multiplication.
\item Least Square Regression
\end{enumerate}

As usual, we will make our lives easier by considering approximate variants (of the optimization versions) of the above problems.


\section{Matrix Multiplicaiton}

\begin{definition}
\emph{(Exact) Matrix Multiplication} is the following problem:
\begin{itemize}
\item Given $A,B \in \mathbb{R}^{n \times d}$,
\item Compute: $C = A^{\top}B \in \mathbb{R}^{d \times d}$.
\end{itemize}
\end{definition}

In general, you may consider the problem for arbitrary fields $\mathcal{K}$,
but we will restrict our attention to $\mathbb{R}$.
(One may also consider matrices of arbitrary dimension.)

Naively, we can solve the above problem in time $O(nd^2)$.
The state of the art for $n \times n$ matrices is time $O(n^\omega)$ for $\omega \approx 2.36\ldots$.
This will yield an algorithm for our problem with complexity $O(d^2n^{\omega-2})$.

However as usual, we are interested in a near linear time, $\sim O(nd)$, algorithm.
To do this exactly is hard, so we will relax the problem to an approximate version.

First, we define the following norm to characterize our approximation guarrantee;
\begin{definition}
  For a matrix $Z \in \mathbb{R}^{m \times n}$, the \emph{(squared) frobenius norm} is defined as follows:
  \[ \|Z\|_F^2 := \sum_{i,j} Z_{i,j}^2. \]
\end{definition}

\begin{definition}
  \emph{(Approximate) Matrix Multiplication} is the following problem:
  \begin{itemize}
  \item Given $A,B \in \mathbb{R}^{n \times d}$,
  \item Compute: $C' \in \mathbb{R}^d$ such that the following holds with high probility,
    \[\|C'- A^{\top}B\|_F \leq \varepsilon \|A\|_F\times\|B\|_F.\]
  \end{itemize}
\end{definition}

Some notation for what follows:
\begin{equation*}
  A = \left[\begin{array}{c} x_1^\top \\ \vdots \\ x_n^\top \end{array}\right]
  \qquad
  B = \left[\begin{array}{c} y_1^\top \\ \vdots \\ y_n^\top \end{array}\right]
\end{equation*}

\subsection{A First Algorithm: Sampling via a Horovitz-Thompson Estimator}

We begin by noting the following:
\begin{claim}
  $A^\top B= \sum_{k=1}^n x_k y_k^\top$ ($xy^\top$ is the ``outer-product'' of vectors $x$ and $y$).
\end{claim}
\begin{proof}
  \[C_{ij} = \left(\sum_{k=1}^n x_k y_k^t\right)_{ij} = \sum_{k=1}^n x_{ki}y_{kj}.\]
\end{proof}

From this, we derive the following algorithm (we will fix parameters in the analysis):
\begin{itemize}
\item Sample $m$ coordinates $k_t$ from $[n]$ 
  ($2m$ vectors: $x_{k_t},y_{k_t}$, $t\in[m]$)
  where the probability of sampling coordinate $k$ is
  $p_k \propto \|x\|_k \|y\|_k$.
\item Then simply output,
  \[ C' = \sum_{t=1}^m \frac{x_{k_t} y_{k_t}}{p_{k_t}}.\]
\end{itemize}

\begin{theorem}
  \[ \Pr\left[ \| C' - C\|_F > \varepsilon \|A\|_F \|B\|_F\right]< \frac{1}{\varepsilon^2 m}.\]
\end{theorem}
Notice that this means we can take $m = \Omega(1/\varepsilon^2)$.

\begin{proof}

\begin{itemize}
\item \emph{Expectation}

  \begin{align*}
    \E[C'] &= \frac{1}{m}\E\left[\sum_{t=1}^m \frac{x_{k_t}y_{k_t}^\top}{p_{k_t}}\right]\\
    &= \frac{1}{m}\sum_{t=1}^m\sum_{k=1}^n \frac{p_k x_k y_k^\top}{p_k}\\
    &= \sum_{k=1}^n x_k y_k^\top = C.
  \end{align*}
  
\item \emph{Variance} %Is this a standard characterization of matrix variance?
  
  \begin{align*}
    V & = \E\left[ \|C' - C\|_F^2\right] \\
    &= \E\left[\sum_{i,j}(C'_{ij}-C_{ij})^2\right]\\
    &= \sum_{i,j} \var[C'_{ij}]\\
    &\leq \sum_{i,j} \var \left[\tfrac{1}{m}\sum_{t=1}^m \underbrace{\frac{x_{k_t i}y_{k_t j}}{p_{k_t}}}_{\mbox{id. dist. var.}}\right]\\
    &= \sum_{i,j} \tfrac{1}{m}\var\left[\frac{x_{ki}y_{kj}}{p_k}\right] \qquad \mbox{(randomness over $k$)}\\
    &\leq \tfrac{1}{m} \sum_{i,j} \E\left[ \left(\frac{x_{ki}y_{kj}}{p_k}\right)^2\right]\\
    &= \tfrac{1}{m} \sum_{i,j} \sum_{k=1}^n p_k \left(\frac{x_{ki}y_{kj}}{p_k}\right)^2\\
    &=\tfrac{1}{m}\sum_{k=1}^n \frac{1}{p_k} \sum_{i,j} x_{ki}^2 y_{kj}^2\\
    &=\tfrac{1}{m} \sum_{k=1}^n \frac{1}{p_k} \|x_k\|_F^2\|y_k\|_F^2
  \end{align*}

  So, take
  \[ p_k := \frac{\|x_k\|_F\|y_k\|_F}{\sum_{i=1}^n \|x_i\|_F \|y_i\|_F}.\]

  Then (via Cauchy-Schwartz),
  \begin{equation*}
    V \leq \frac{\left(\sum_{k=1}^n \|x_k\|_F\|y_k\|_F\right)^2}{m}
    \leq \frac{\left(\sum_{k=1}^n \|x_k\|_F^2 \right)\left(\sum_{k=1}^n \|y_k\|_F^2\right)}{m}
    = \frac{\|A\|_F^2\|B\|_F^2}{m}
  \end{equation*}
  
\item So applying Chebyshev to the above,

  \begin{equation*}
    \Pr\left[ \|C'-C\|_F^2 > \varepsilon^2 \|A\|_F^2\|B\|_F^2\right] 
    \leq \frac{\E\left[\|C'-C\|_F\right]}{\varepsilon^2\|A\|_F^2\|B\|_F^2} 
    \leq \frac{1}{m\varepsilon^2}
  \end{equation*}
\end{itemize}
\end{proof}

\subsection{A Second Algorithm: Using Dimension Reduction}

\begin{note}
  We can view the above algorithm as the following:
  \begin{itemize}
  \item Choose a random $\Pi \in \mathbb{R}^{m \times n}$ where
    \[ \Pi_{i,j} := \left\{ \begin{array}{rl} \frac{1}{\sqrt{mp_k}} & \mbox{if } (i,j) = (t,k_t)\\ 0  & \mbox{otherwise}\end{array}\right.\]
  \item Compute:
    \[ C' = (\Pi A)^\top(\Pi B).\]
  \end{itemize}
  
  Observe that the above algorithm requires two passes over the data, one to sample $\Pi$ (compute the $p_k$'s)
  and one to compute the ``reduced'' matrix product (or the sum in our previous formulation).
\end{note}

Given this ``randomized-projection/embedding'' formulation of our approximation algorithm,
it seems an appropriate place to invoke the magic of Johnson-Lindenstrauss.
Consider the following definition:

\begin{definition}
  $\Pi \in \mathbb{R}^{m \times n}$ is an \emph{$(\varepsilon,\delta)$-dimension reducing matrix}, $(\varepsilon,\delta)$-DR,
  if
  \[ \forall x\in\mathbb{R}^n, \Pr\left[|\|\Pi x\|_2^2-\|x\|_2^2| > \varepsilon \|x\|_2^2\right]\leq \delta.\]
\end{definition}


Given some $(\varepsilon,\delta)$-DR matrix $\Pi$, our algorithm is to simply compute:
\[ C' = \left( \Pi A\right)^\top \left( \Pi B \right).\]

\begin{theorem}
  $\Pi$ is $(\varepsilon,\delta)$-DR $\implies \Pr\left[\|C'-C\|_F > 3\varepsilon\|A\|_F\|B\|_F\right]\leq 3d^2\delta$.
\end{theorem}

\begin{remark}
  With a more precise version of the Johnson-Lindenstrauss lemma we can remove the $d^2$ factor from the above.
\end{remark}

\begin{corollary}
  If we choose $m=O(1/\varepsilon^2 \log(1/\delta)), \delta = \frac{1}{10d^2}$,
  then (naively) we can compute $C'$ in time $O(mnd)+O(dmd) = O(\frac{nd+d^2}{\varepsilon^2}\log d)$.
  
  By the above remark, the $\log d$ factor is simply an artifact of our analysis.
\end{corollary}

To prove the theorem, we will show $C'_{ij} \approx C_{ij}$ with probability $\geq 1-3\delta$ and then take a union bound (hence the $d^2$).

\begin{proof}
  First some notation:
  \begin{equation*}
    A = \left[\begin{array}{ccc} A_1 & \cdots & A_d \end{array}\right]
    \qquad
    B = \left[\begin{array}{ccc} B_1 & \cdots & B_d \end{array}\right]
  \end{equation*}
  
  \begin{equation*}
    a_i := \frac{A_i}{\|A_i\|_2} \qquad b_i := \frac{B_i}{\|B_i\|_2}
  \end{equation*}

  Note:
  \begin{itemize}
  \item $C_{ij} = A_i^\top B_j = \|A_i\|\|B_j\|a_i^\top b_j$.
  \item With probability $\geq 1-3\delta$,
    \begin{align*}
      C'_{ij} = \left(\Pi A_i\right)^\top \left(\Pi B_j \right) &= \|A_i\|\|B_j\|(\Pi a_i)^\top (\Pi b_j)\\
      &= \|A_i\|\|B_j\| \left[\|\Pi a_i\|^2 + \|\Pi b_j\|^2 - \tfrac{1}{2}\|\Pi a_i - \Pi b_J\|^2\right]\\
      &= \|A_i\|\|B_j\| \left[\|a_i\|^2 +\|b_j\|^2 - \|a_i-b_j\|^2 \pm 3\varepsilon\right] \qquad (\Pi \mbox{ is }(\varepsilon,\delta)-\mbox{DR})\\
      &= \|A_i\|\|B_j\| \left[a_ib_j \pm 3\varepsilon\right]
    \end{align*}
  \end{itemize}
  
  So with probability $\geq 1-3\delta$, $(C'_{ij}-C_{ij})^2 \leq \|A_i\|_2^2\|B_j\|_2^2(3\varepsilon)^2$.
  This implies (via union bound) that with probability $\geq 1-3\delta d^2$,
  \[ \|C'-C\|_F \leq \sum_{ij}9\varepsilon^2\|A_i\|^2\|B_j\|^2 = 9\varepsilon^2\|A\|_F^2\|B\|_F^2.\]
\end{proof}

\section{Least Squares Regression}

\begin{definition}
  \emph{(Exact) Least Squares Regression} is the following problem:
  \begin{itemize}
  \item Given $A\in\mathbb{R}^{n\times d}, b\in \mathbb{R}^n$,
  \item find $x^* = \argmin_{x\in\mathbb{R}^d} \|Ax-b\|_2$.
  \end{itemize}
\end{definition}

We can consider least squares regression as a simple learning problem where
the $i$-th row of $A$, $a^{(i)}$, is labeled with $b_i$ according to some approximately linear function.
\begin{definition}
  A function $f:\mathbb{R}^d \to \mathbb{R}$ is \emph{linear} if
  \[\exists y\in\mathbb{R}^d: f(x)= \langle x,y \rangle.\]
\end{definition}

If $\exists x: Ax=b$ then the problem is easy.
In general, we are only assume $\exists x: Ax\approx b$.

In general, we can do least squares regression via Singular Value Decomposition (in time $\tilde{O}(nd^{\omega-1})$),
but perhaps we can speed things up by loosening the approximation.

\begin{definition}
  \emph{(Approximate) Least Squares Regression} is the following problem:
  \begin{itemize}
  \item Given $A\in\mathbb{R}^{n\times d}, b\in \mathbb{R}^n$,
  \item Let $x^* = \argmin_{x\in\mathbb{R}^d} \|Ax-b\|_2$.
    Find $x\in\mathbb{R}^d$ such that
    \[ \|Ax-b\|_2 \leq (1+\varepsilon)\|Ax^*-b\|_2.\]
  \end{itemize}
\end{definition}

To solve this problem we will use dimension reduction, as promised.

First, we define a special kind of dimension reducing matrix:
\begin{definition}
  $\Pi \in \mathbb{R}^{m \times n}$ is a \emph{$(d,\varepsilon,\delta)$-subspace embedding}, 
  $(d,\varepsilon,\delta)$-SE,
  if $\forall P \subset \mathbb{R}^n$ such that $P$ is a $d$-dimensional subspace,
  \[\Pr[\forall p \in P: |\|\Pi p\| - \|p\||\leq \varepsilon \|p\|]\geq 1-\delta.\]
\end{definition}

Then given some such SE $\Pi$, our alorithm is simply:
find $\argmin_{x}\| \Pi A x - \Pi b \|$ (via SVD).

Naively,
the time to reduce dimension is $O(mnd)$.
The time to perform SVD on the result is $O(md^{\omega-1}$.

So if we take $m=O(d/\varepsilon^2)$,
then the resulting algorithm has time complexity
\[ O(\frac{nd^2}{\varepsilon}+md^{\omega-1}).\]

If we use a faster version of Johnson-Lindenstrauss, 
we can acheive $O_\varepsilon\left((n\log n + m^3)d\right)$ time complexity.

Unfortunately, at this point we ran out of time.
We will finish up this application of dimension reduction next lecture.



\iffalse
\begin{thebibliography}{9}
\end{thebibliography}
\fi

\end{document}