Streaming PCA: Matching Matrix Bernstein and Near-Optimal Finite Sample Guarantees for Oja's Algorithm
Prateek Jain\footnote{Microsoft Research India. Email: prajain@microsoft.com
This work provides improved guarantees for % solving a key problem in data analysis, streaming principle component analysis (PCA). Given $\A_1, \ldots, \A_n\in \R^{d\times d$ sampled independently from distributions satisfying $\E{\A_i = \Cov$ for $\Cov \succeq \mZero$, this work provides an $O(d)$-space linear-time single-pass streaming algorithm for estimating the top eigenvector of $\Cov$. The algorithm nearly matches (and in certain cases improves upon) the accuracy obtained by the standard batch method that computes top eigenvector of the empirical covariance $\frac{1{n \sum_{i \in [n] \A_i$ as analyzed by the matrix Bernstein inequality. Moreover, to achieve constant accuracy, our algorithm improves upon the best previous known sample complexities of streaming algorithms by either a multiplicative factor of $O(d)$ or $1/\mathrm{gap$ where $\mathrm{gap$ is the relative distance between the top two eigenvalues of $\Cov$. These results are achieved through a novel analysis of the classic Oja's algorithm, one of the oldest and most popular algorithms for streaming PCA. In particular, this work shows that simply picking a random initial point $\w_0$ and applying the update rule $\w_{i + 1 = \w_i + \eta_i \A_i \w_i$ suffices to accurately estimate the top eigenvector, with a suitable choice of $\eta_i$. We believe our result sheds light on how to efficiently perform streaming PCA both in theory and in practice, and we hope that our analysis may serve as the basis for analyzing many variants and extensions of streaming PCA.
approach.tex
\section{Approach}\label{sec:approach}

Let us now describe the approach to analyze Oja's algorithm. 
We provide our main theorem regarding the convergence rate of Oja's algorithm and discuss how it is proved. 
The details of the proof are deferred to Section~\ref{sec:proofs} and the use of the theorem to choose step sizes is in Section~\ref{sec:results}.

One of the primary difficulties in analyzing Oja's algorithm, or more broadly any algorithm for streaming PCA, is choosing a subtle potential function to analyze the  method. If we try to analyze the progress of Oja's algorithm in every iteration $i$, by measuring the quality of $\w_i$, we run the risk that during the first few iterations of Oja's algorithm a step may actually yield a $\w_{i + 1}$ that is orthogonal to $\v_{i}$. If this happens, even in the typical best case, where all future samples are $\Cov$ itself, we would still fail to converge. In short, if we do not account for the randomness of $\w_0$ in our potential function then it is difficult to show that a rapidly convergent algorithm does not catastrophically fail.

%To overcome this difficulty, 
Rather than analyzing the convergence of $\w_i$ directly we instead analyze the convergence of Oja's algorithm as an operator on $\w_0$.  Oja's algorithm simply considers the matrix 
\begin{equation}
\label{eq:onestep-1}
\B_{n} \defeq (\eye + \eta_n \A_n )(\eye + \eta_{n - 1} \A_{n - 1}) \cdots (\eye + \eta_1 \A_1)
\end{equation}
and outputs the normalized result of applying this matrix, $\B_{n}$, to the random initial vector, i.e.
\begin{equation}
\label{eq:onestep-2}
\w_n  = \frac{\B_n \w_0}{\| \B_n \w_0 \|_2} ~.
\end{equation}
Rather than analyze the improvement of $\w_{n + 1}$ over $\w_n$ we analyze $\B_{n + 1}$'s improvement over $\B_n$.

Another interpretation of \eqref{eq:onestep-1} and \eqref{eq:onestep-2} is that Oja's algorithm simply approximates $\v_n$ by performing 1 step of the power method on the matrix $\B_{n}$. Fortunately, analyzing when 1 step of the power method succeeds is fairly straightforward as we show below:

\begin{lemma}[One Step Power Method]\label{lem:onesteppowermethod}
Let $\B \in \R^{d \times d}$, let $\vtilde \in \R^{d}$ be a unit vector, and let $\Vtildeperp$ be a matrix whose columns form an orthonormal basis of the subspace orthogonal to $\vtilde$. If $\w \in \R^{d}$ is chosen uniformly at random from the surface of the unit sphere then with probability at least $1-\delta$
	\begin{align*}
	\sin^2\left(\vtilde, \frac{\B\w}{\|\B\w\|_2}\right)=1-\left( \frac{\vtilde^\top\B \w}{\twonorm{\B \w}}\right)^2\leq \frac{C\log \left(1/\delta\right)}{\delta}\frac{\trace{\VtildeperpTr\B\trans{\B}\Vtildeperp}}{\trans{\vtilde}\B\trans{\B}\vtilde}
	\end{align*}
where $C$ is an absolute constant.
\end{lemma}

\begin{proof}
%Let $\B=\V_{\B} \mSigma_{\B} \U_{\B}^\top $ be the (full)singular value decomposition of $\B$. Also, let $\valpha=\U_{\B}^\top\w$, i.e., $w=\U_{\B}\valpha$. Note that $\valpha$ is also distributed uniformly over the sphere, i.e., $\valpha=\vg/\|\vg\|_2$ where $\vg\sim N(0, I)$.  
As $\w$ is distributed uniformly over the sphere, we have: $\w=\vg/\|\vg\|_2$ where $\vg\sim N(0, I)$. Consequently, with probability at least $1-\delta$
\begin{align*}
1-\left(\frac{\vtilde^\top\B \w}{\twonorm{\B \w}}\right)^2&=\frac{\vg^\top \B^\top (\eye-\vtilde\vtilde^\top)\B\vg}{\vg^\top \B^\top \B \vg}\stackrel{\zeta_1}{\leq}\frac{C_1}{\delta} \frac{\vg^\top \B^\top (\eye-\vtilde\vtilde^\top)\B\vg}{\vtilde^\top\B\B^\top\vtilde}\\
&\stackrel{\zeta_2} {\leq}\frac{C\log(1/\delta)}{\delta} \frac{\trace{\B^\top (\eye-\vtilde\vtilde^\top)\B}}{\vtilde^\top\B\B^\top\vtilde},
\end{align*}
where $C_1$ and $C$ are absolute constants. $\zeta_1$ follows as $\vg^\top\B^\top\B\vg\geq (\vtilde^\top\B\vg)^2\geq \frac{1}{C_1}\vtilde^\top\B\B^\top\vtilde$ where the second inequality follows from the fact that $\vtilde^\top\B\vg$ is a Gaussian random variable with variance $\twonorm{\B^\top \vtilde}^2$. Similarly, $\zeta_2$ follows from the fact that$\vg^\top \B^\top (\eye-\vtilde\vtilde^\top)\B\vg$ is a $\chi^2$ random variable with $\trace{\B^\top (\eye-\vtilde\vtilde^\top)\B}$-degrees of freedom. 
% \[
% \frac{\B w}{\twonorm{\B w}}
% = \frac
% {\sum_{i \in [n]} \alpha_i \sigma_{\B}(i) \v_{\B}(i)}
% {\sqrt{\sum_{i \in [n]} \alpha_i^2 \sigma_{\B}(i)^2}}
% \]
% Where $\sigma_{\B}(i) = [\mSigma_{\B}]_{ii}$ and $\v_{\B}(i)$ is row $i$ of $\V_{\B}$. Consequently,
% \[
% \left\langle
% \utilde , \frac{\B w}{\twonorm{\B w}}
% \right\rangle^2
% \]
% ... this is enough to see why true. TODO Should write out rest of proof  once decide how want to phrase statement
\end{proof}

This lemma makes our goal clear. To show that Oja's algorithm succeeds we simply need to show that with constant probability $\v_1^\top \B_n \B_n^\top \v_1$ is relatively large and $\trace{\Vperp \B_n \B_n^\top \Vperp}$ is relatively small, where $\Vperp$ is a matrix whose columns form an orthonormal basis of the subspace orthogonal to $\v_1$. This immediately alleviates the issues of catastrophic failure that plagued analyzing $\w_n$. So long as we pick $\eta_i$ sufficiently small, i.e. $\eta_i = O(1/\max\{\upperbound,\lambda_1\})$ then $\eye + \eta_i \A_i$ is invertible. In this case $\B_n \B_n^\top$ is invertible and $\v_1^\top \B_n \B_n^\top \v_1 > 0$. In short, so long as we pick $\eta_i$ sufficiently small the quantity we wish to bound $\trace{\Vperp \B_n \B_n^\top \Vperp} / \v_1^\top \B_n \B_n^\top \v_1$ is always finite.

To actually bound $\v_1^\top \B_n \B_n^\top \v_1$ and
$\trace{\Vperp^{\top} \B_n \B_n^\top \Vperp}$ we split the analysis
into several parts in Section~\ref{sec:proofs}. First, we show that
$\E{\trace{\Vperp^{\top} \B_n \B_n^\top \Vperp}}$ is small, which
implies by Markov's inequality that $\trace{\Vperp^{\top} \B_n
  \B_n^\top \Vperp}$ is small with constant probability. Then, we show
that $\Eop \v_1^\top \B_n \B_n^\top \v_1$ is large and that $\Var{
  \v_1^\top \B_n \B_n^\top \v_1}$ is small. By Chebyshev's inequality
this implies that ${\v_1^\top \B_n \B_n^\top \v_1}$ is large with
constant probability. Putting these together we achieve the main
technical result regarding the analysis of Oja's method. Once we
devise this roadmap, the proof is fairly straightforward.  

\begin{theorem}[Oja's Algorithm Convergence Rate] \label{thm:main1}
Let $\delta > 0$ and step sizes $\eta_i \leq \frac{1}{4\cdot \max\{M,\lambda_1\}}$. 	The output $\wt{n}$ of Algorithm~\ref{algo:oja} is an $\epsilon$-approximation to $\v_1$ with probability at least  $1 - \delta$ where
\begin{align*}
\epsilon &\leq \frac{1}{Q} \exp\left(5 \overline{\variance}\sum_{i\in[n]}\eta_i^2\right)
\left(d \cdot \exp\left(-2(\lambda_1-\lambda_2)\sum_{i \in [n]} \eta_i  \right) + \variance \sum_{i=1}^{n} \eta_i^2 \exp
\left(-\sum_{j =i+1 }^n 2\eta_j(\lambda_1-\lambda_2) \right)\right)
\end{align*}
 where $Q \eqdef \frac{\delta^2}{C\log(1/\delta)}\left( 1 - \frac{1}{\sqrt{\delta}}\sqrt{ \exp\left(18\overline{\variance}\sum_{i=1}^{n} \eta_i^2 \right) -  1}\right)$, $\overline{\variance} \eqdef \variance+\lambda_1^2$, and $C$ is an absolute constant.
\end{theorem}

Theorem~\ref{thm:main1} is proved in Section~\ref{sec:proofs}. Theorem~\ref{thm:main1} serves as the basis for our results regarding Oja's algorithm. In the next section we show how to use this theorem to choose step sizes and achieve the main results of this paper. 
conclusion.tex
\section{Conclusion and Future Work} \label{sec:conc}
This work presented a finite sample complexity and asymptotic convergence rates for the classic Oja's algorithm for top-$1$ component streaming PCA that match well known matrix concentration and perturbation results for computing the top eigenvector. In fact, asymptotically our bound improves upon standard matrix Bernstein bounds by a factor of $\order{\log d}$. Our results are tighter than existing streaming PCA results by a factor of either $\order{d}$ or $\order{1/\textrm{gap}}$.%; getting rid of $\log d$ in matrix Bernstein is a well-known open problem. 

Our analysis relied on a novel view of the algorithm and is technically fairly simple. We hope that our analysis opens a way to make progress on the many variants of PCA that occur in both theory and practice.
%such as results free of eigenvalue gap, optimal success probability, extension to more components and so on. We hope that the simplicity of our analysis promotes further research into this important problem.
In particular, we believe the following directions should be of wide interest:
\begin{itemize}
	\item \textbf{Multiple components}: Currently, our result holds only for estimating the top eigenvector of $\Cov$. Extension of our technique to compute top-$k$ eigenvectors is an important future direction. 
	\item \textbf{Rayleigh quotient}: Another standard metric to measure optimality of $\wt{n}$ is Rayleigh quotient: $\wt{n}^\top \Cov \wt{n}$. Converting our bounds on $\sin^2(\wt{n},\v_1)$ to Rayleigh quotient loses a multiplicative factor of $\order{1/\textrm{gap}}$ compared to the optimal rate. A direct analysis that does not lose this factor is an interesting open problem. Results on Rayleigh quotient may also help in obtaining sample complexity guarantees that are independent of eigenvalue gap. %	\item \textbf{Gap Free Bounds}: Can we get error guarantees even without the assumption of an eigenvalue gap, i.e. $\lambda_1 > \lambda_2$, e.g. in terms of Rayleigh Quotient error.
	\item \textbf{High Probability}: This work focused on obtaining tight bounds on the error. However, the dependence of our results on success probability is quite suboptimal. One way to fix this is to run many copies of the algorithm, each with say $3/4$ success probability and then output the geometric median of the solutions, which can be done in nearly linear time\cite{cohen2016median}. However, we conjecture that a tighter analysis using our techniques might directly lead to improved dependency on success probability and possibly help solve some of the other problems mentioned above. %Can we succeed with high probability (without resorting to median trick, i.e. make it more natural). Also can we make the constant success only due to the random initial, not the operator.
\end{itemize}

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

intro.tex

\section{Introduction}\label{sec:intro}

Principal component analysis (PCA) is one of the most fundamental problems in machine learning, numerical linear algebra, and data analysis. It is commonly used for data compression, image processing, and visualization  \cite{jolliffe2002principal} etc.

When we desire to perform PCA on large data sets, it may be the case that we cannot afford more than single pass over the data (or worse to even store the data in the first place) \cite{hall1998incremental,weng2003candid,ross2008incremental}. To alleviate this issue, a popular line of research over the past several decades has been to consider streaming algorithms for PCA under the assumption that the data has reasonable statistical properties \cite{krasulina1970method,oja1982simplified,BalsubramaniDF13,MitgliakisCJ13,SaOR15}. There have been significant breakthroughs in getting near-optimal streaming PCA algorithms under fairly specialized models, e.g. spiked covariance \cite{SaOR15}. 

This work considers one of the most natural variants of PCA,
estimating the top eigenvector of a symmetric matrix, under a mild
(and standard) set
of assumptions under which concentration of measure applies (under the
matrix Bernstein
inequality\cite{vershynin2010introduction,Tropp12}). In particular,
the setting is as follows: 

\begin{definition}[Streaming PCA]
\label{def:streaming-PCA}
Let $\A_1, \A_2, ... , \A_n \in \R^{d \times d}$ be a sequence  of (not necessarily symmetric) matrices sampled independently from distributions that satisfy the following: \begin{enumerate}
	\item	$\E{\At{i}} = \Cov$ for symmetric positive semidefinite (PSD) matrix $\Cov \in \R^{d \times d}$,
	\item	$\twonorm{\At{i} - \Cov} \leq \upperbound$ with probability $1$, and
	\item	$\max\left\{\twonorm{\E{(\At{i}-\Cov)\trans{(\At{i} - \Cov)}}}, \twonorm{\E{\trans{(\At{i}-\Cov)}(\At{i}-\Cov)}}\right\} \leq \variance$.
\end{enumerate}
Let $\v_1, ..., \v_d$ denote the eigenvectors of $\Cov$ and $\lambda_1 \geq ... \geq \lambda_d$ denote the corresponding eigenvalues. Our goal is to compute an $\epsilon$-approximation to $\v_1$, that is a unit vector $\w$ such that
$
\sin^2(\w, \v_1) \eqdef 1 - (\w^\top \v_1)^2 \leq \epsilon
$, in a single pass while minimizing space, time, and error (i.e. $\epsilon$). Note that $\sin(\w, \v_1)$ denotes the $\sin$ of the angle between $\w$ and $\v_1$.
\end{definition}

A special case of Streaming PCA is to estimate the top eigenvector of the covariance matrix of a distribution $\dist$ over $\R^{d}$, i.e.  given independent samples $\va_1, ..., \va_n \in \R^d$ estimate the top eigenvector of $\Eop_{\va \sim \dist}[\va \va^\top]$. This encompasses the popular "spiked covariance model"\cite{johnstone2001distribution}.

It is well known that to solve the Streaming PCA problem, one can simply compute the empirical covariance matrix $\frac{1}{n} \sum_{i \in [n]} \A_i$ and compute the right singular vector of this matrix. Here, matrix Bernstein inequality\cite{vershynin2010introduction,Tropp12} and Wedin's theorem\cite{wedin1972perturbation} implies the following standard sample complexity bound for the Streaming PCA problem:

\begin{theorem} [Eigenvector Concentration using matrix Bernstein and Wedin's theorem]
\label{thm:bernstein_kahan}
Under the assumptions of Definition~\ref{def:streaming-PCA}, 
%\begin{equation*}
% \ge \max\left\{ \frac{16 \variance}{(\lambda_1 - \lambda_2)^2}\log \frac{d}{\delta}, 
%\frac{4\upperbound}{(\lambda_1 - \lambda_2) }\log\frac{d}{\delta}\right\}.
%\end{equation*}
the top right singular vector $\widehat{\v}$ of $\widehat{\Cov}=\frac{1}{n} \sum_{i \in [n]} \A_i$ is an $\epsilon$-approximation to the top eigenvector $\v_1$ of $\Cov$ with probability $1-\delta$, where  % for $\epsilon$ satisfying
\begin{equation*}
\sin^2(\widehat{\v}, \v_1) \leq \epsilon
\leq \frac{16 \variance \log \frac{d}{\delta}}{(\lambda_1 -\lambda_2)^2} \cdot \frac{1}{n} + \left(\frac{4\upperbound \log\frac{d}{\delta}}{\lambda_1 - \lambda_2 }\right)^2 \cdot \frac{1}{n^2}.
\end{equation*}
\end{theorem}

Theorem~\ref{thm:bernstein_kahan} is essentially the previous best
sample complexity known for estimating the top eigenvector~\footnote{In recent work in\cite{jin2015robust} it was shown
  that the $log(d/\delta)$ factor in the first term could be removed
  asymptotically for small enough $\epsilon$ if only constant success
  probability is required.}. 
%for estimating the top eigenvector 
Unfortunately, the above is purely a statistical claim, and,
algorithmically, there are least two concerns. First, computing the
empirical covariance matrix
$\widehat{\Cov}=\frac{1}{n}\sum_{i \in [n]} \A_i$ naively requires
$O(d^2)$ time and space, and second, computing the top eigenvector of
the empirical covariance matrix in general may require super linear
time\cite{golub2012matrix}. While there have been many attempts to
produce streaming algorithms that use only $O(d)$ space to solve the
streaming PCA problem, to our knowledge, all previous methods
either lose a multiplicative factor of either
$\frac{\lambda_1}{\lambda_1 - \lambda_2}$ or $d$ in the analysis in
order to achieve constant accuracy when applied in our
setting\cite{BalsubramaniDF13,MitgliakisCJ13, HardtP14,
  SaOR15,jin2015robust}.

In an attempt to overcome this limitation and improve the guarantees
for solving the streaming PCA problem, this work seeks to address the following question:

\begin{quote}
\emph{
Can we match the sample complexity of matrix Bernstein + Wedin's theorem with an algorithm that uses $O(d)$ space only and takes a single linear-time pass over the input?}
\end{quote}

This work answers this question in the affirmative, showing that one
can succeed with constant probability matching the sample complexity
of Theorem~\ref{thm:bernstein_kahan} up to logarithmic terms and small
additive factors. Interestingly, this is achieved by providing a novel
analysis of the classical Oja's algorithm, which is perhaps, the most
popular algorithm for Streaming PCA\cite{oja1982simplified}.

\begin{algorithm}[H]
	\caption{Oja's algorithm for computing top eigenvector}
	\begin{algorithmic}\label{algo:oja}
		\renewcommand{\algorithmicrequire}{\textbf{Input: }}
		\renewcommand{\algorithmicensure}{\textbf{Output: }}
		\REQUIRE $\At{1},\cdots,\At{n}$.
		\STATE Choose $\wt{0}$ uniformly at random from the unit sphere
		\FOR{$t = 1,\cdots,n$}
		\STATE $\wt{i} \leftarrow \wt{i-1} + \eta_i \At{i}\wt{i-1}$
		\STATE $\wt{i} \leftarrow \wt{i}/\twonorm{\wt{i}}$
		\ENDFOR
		\ENSURE  $\wt{n}$
%		\STATE \textbf{Return} $\wt{n}$.
	\end{algorithmic}
\end{algorithm}

Oja's algorithm is one of the simplest algorithms one would imagine
for the streaming PCA problem (See Algorithm~\ref{algo:oja}). In fact,
due to its simplicity, it was proposed a neurally plausible algorithm.
In the case that each $\A_i$ comes from the same distribution $\dist$
it corresponds to simply performing projected stochastic gradient
descent on the objective function of maximizing the Rayleigh Quotient
over the distribution
$\max_{\twonorm{w} = 1} \Eop_{\A \sim \dist} \w^\top \A \w$. It is
well known that under very mild conditions on the stepsize sequence,
Oja's algorithm asymptotically converges to the top eigenvector of the
covariance matrix $\Cov$\cite{oja1982simplified}. However, obtaining
optimal rates of convergence, let alone finite sample guarantees, for
Streaming PCA has been quite challenging. The best known results are
off from Theorem~\ref{thm:bernstein_kahan} by a factor of
$\order{d}$\cite{SaOR15}.% or $\order{\frac{\lambda_1}{\lambda_1-\lambda_2}}$\cite{HardtP14}.

This work shows that for proper choice of learning rates $\eta_i$,
Oja's algorithm in fact can improve the best known results for
streaming PCA and answer our question in the affirmative. In
particular, we have that:

\begin{theorem}\label{thm:intro1}
Let the assumptions of Definition~\ref{def:streaming-PCA} hold. Suppose the step size sequence for Algorithm~\ref{algo:oja} is chosen to be $\eta_i = \frac{\log d}{(\lambda_1-\lambda_2)(\beta+i)}$, where
	\begin{align*}
	\beta \eqdef 40\max\left(\frac{ \upperbound \log d}{(\lambda_1-\lambda_2)},\frac{\left(\variance+\left(\lambda_{1}\right)^2\right) \log^2 d}{(\lambda_1-\lambda_2)^2}\right).
	\end{align*}
	Then the output $\wt{n}$ of Algorithm~\ref{algo:oja} is an $\epsilon$-approximation to the top eigenvector $\v_1$ of $\Cov$ satisfying
	\begin{align*}
	\sin^2(\w_n,\v_1)\leq \epsilon \;\leq \;C \left(\frac{\variance \log d}{(\lambda_1-\lambda_2)^2}\cdot \frac{1}{n} + \left(\frac{2\beta}{n}\right)^{2 \log d}\right),
	\end{align*}
	with probability greater than $3/4$. Here $C$ is an absolute numerical constant.
\end{theorem}
The error above should be interpreted as being the sum of a
$\order{\frac{1}{n}}$ higher order term and another
$\order{(2\beta/n)^{2 \log d}}$ lower order term which is at most
$o(\frac{1}{n^{\log d}})$ (once $n > 4\beta^2$). In particular, this
result shows that, up to an additive lower order term, one can match
Theorem~\ref{thm:bernstein_kahan} with an asymptotic error of
$\order{\frac{\variance \log d}{(\lambda_1-\lambda_2)^2 n}}$ with
constant probability. The lower order term has $\beta$ which is the
$\max$ of three parts:
$\frac{ \upperbound \log d}{(\lambda_1-\lambda_2)}$,
$\frac{ \variance \log^2 d}{(\lambda_1-\lambda_2)^2}$ and
$\frac{ \lambda_1^2 \log^2 d}{(\lambda_1-\lambda_2)^2}$. The first
part, depending on $\upperbound$, is exactly the same as what appears
in Theorem~\ref{thm:bernstein_kahan}. The second one, depending on
$\variance$ has an additional $\log d$ factor over the first order
term and is irrelevant once, say $n > 10 \beta$. Notably, the third part,
depending on $\lambda_1^2$, does not appear in
Theorem~\ref{thm:bernstein_kahan}; it arises here entirely due to
computational reasons: the setting allows only a \emph{single linear-time
  pass} over the matrices, while Theorem~\ref{thm:bernstein_kahan}
makes no such assumption. For instance, consider the case
$\variance=0$ which means $\At{1}=\Cov$. Matrix Bernstein tells us
that one sample is sufficient to compute $\v_1$. However, it is not
evident how to compute it using a single pass over $\At{1}$. Note
however, that the rate at which the lower order terms,
i.e. $o\left(\frac{1}{n^{\log d}}\right)$, decrease is much better than
$\order{1/n^2}$ guaranteed by Theorem~\ref{thm:bernstein_kahan}.

In fact, this result also improves the asymptotic error rate obtained
by Theorem~\ref{thm:bernstein_kahan}. In particular, the
following result shows that Oja's algorithm gets an asymptotic rate of
$\order{\frac{\variance }{(\lambda_1-\lambda_2)^2 n}}$ which is better
than that of matrix Bernstein by a factor of
$\order{\log d}$.\footnote{A similar asymptotic result was recently
  obtained by\cite{jin2015robust}. However, their result requires an
  initial vector that is constant close to $\v_1$, which itself is a
  difficult problem.} \begin{theorem}\label{thm:intro2} Let the
  assumptions of Definition~\ref{def:streaming-PCA} hold. Suppose the
  step size sequence for Algorithm~\ref{algo:oja} is chosen to be
  $\eta_i = \frac{6}{(\lambda_1-\lambda_2)(\beta+i)}$, where
  \begin{align*} \beta \eqdef 720\max\left(\frac{ \upperbound
      }{(\lambda_1-\lambda_2)},\frac{\variance + \lambda_{1}^2
      }{(\lambda_1-\lambda_2)^2}\right). \end{align*} Suppose
  $n > \beta^{1.2}d^{0.1}$. Then the output $\wt{n}$ of
  Algorithm~\ref{algo:oja} is an $\epsilon$-approximation to the top
  eigenvector $\v_1$ of $\Cov$ satisfying \begin{align*}
    \sin^2(\w_n,\v_1)\leq \epsilon \;\leq \; C \left(\frac{ \variance
      }{(\lambda_1-\lambda_2)^2}\cdot \frac{1}{n} +
      \frac{1}{n^2}\right), \end{align*} with probability greater than
  $3/4$. Here $C$ is an absolute numerical constant. \end{theorem}

Note that Theorems~\ref{thm:intro1} and~\ref{thm:intro2} guarantee
success probability of $3/4$. One way to boost the probability to
$1-\delta$, for some $\delta > 0$, is to run $\order{\log 1/\delta}$
copies of the algorithm, each with $3/4$ success probability and then
output the geometric median of the solutions, which can be done in
nearly linear time\cite{cohen2016median}. The detailes are omitted
here.
%Alternately, one This can indeed be improved to $1-\delta$ for any $\delta > 0$, at the expense of a worse constant $C = \order{\log(1/\delta)/\delta^2}$. Please see Theorem~\ref{thm:main} for the full statement. Note that this dependence on success probability is suboptimal compared to that of Theorem~\ref{thm:bernstein_kahan} and closing this gap (or showing that it is not possible) is an interesting direction of future work.

Beyond the improved sample complexities we believe our analysis sheds
light on the type of step sizes for which Oja's algorithm converges
quickly and therefore illuminates how to efficiently perform streaming
PCA. We note that we have essentially assumed an oracle which sets the
step size sequence, and an important question is how to set the step
size in a robust and data data driven manner. Moreover, we believe that our
analysis is fairly general and hope that it may be extended to make
progress on analyzing the many variants of PCA that occur in both
theory and in practice.

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

lemmas.tex
%!TEX root = Ojasalgo.tex


\section{Bounding the Convergence of Oja's Algorithm}
\label{sec:proofs}

In this section, we present a detailed proof of Theorem~\ref{thm:main1}. The proof follows the approach outlined in Section~\ref{sec:approach} and uses the notation of that section, i.e.
\begin{itemize}
	\item We let $\Bt{n} \eqdef \left(\eye + \eta_n \At{n}\right)\cdots \left(\eye + \eta_1 \At{1}\right)$ with $\Bt{0} \eqdef \eye$
	\item We let $\variancetotal \eqdef \variance+\lambda_1^2$
	\item We let $\Vperp \in \R^{d \times d - 1}$ denote a matrix whose columns form an orthonormal basis for the subspace orthogonal to $\v_1$ ~.
\end{itemize}
We first provide several technical lemmas bounding the expected behavior of $\B_n$ and ultimately use these lemmas to prove Theorem~\ref{thm:main1}. We begin with a straightforward lemma bounding the rate of increase of $\E{\B_t \B_t^\top}$ in spectral norm.

\begin{lemma}\label{lem:spectralboundE}
For all $t \geq 0$ and $\eta_i \geq 0$ we have
\[
\twonorm{\E{\Bt{t} \BtTr{t} }} \leq
\exp\left(\sum_{i\in[t]} 2\eta_i \lambda_1 + \eta_i^2 \variancetotal\right)
~.
\]
\end{lemma}

\begin{proof}
Let $\alpha_t \eqdef \|\E{\Bt{t} \Bt{t}^\top} \|_2$, i.e., $\E{\Bt{t} \Bt{t}^\top} \preceq \alpha_t \eye$. For all $t > 0$, 
\begin{align}
\E{\Bt{t} \Bt{t}^\top} 
&=
\E{(\eye + \eta_{t} \A_{t}) \Bt{t - 1} \Bt{t - 1}^\top (\eye + \eta_{t} \A_{t})^\top}
\preceq \alpha_{t-1} \E{ (\eye + \eta_{t} \A_{t}) (\eye + \eta_{t} \A_{t}^\top)},\notag
\\
&= \alpha_{t - 1} \E{\eye + \eta_t \A_t + \eta_t \A_t^\top + \eta_t^2 \A_t \A_t^\top}\preceq \alpha_{t - 1} \left[\eye + 2\eta_t \Cov + \eta_t^2 (\Cov^2 + V \eye) \right],\label{eq:sbE_1}
\end{align} 
where the last inequality follows from $\E{A_t}=\Cov$ and, 
\[
\E{\A_t\A_t^\top} = \Cov^2 + \E{(\A_t - \Cov)(\A_t - \Cov)^\top} \preceq \Cov^2 + V \eye
~.
\]
Using \eqref{eq:sbE_1} along with $\|\E{\Bt{t} \Bt{t}^\top}\|_2=\alpha_t$, $\Cov \preceq \lambda_1 \eye$, and $\Cov^2 \preceq \lambda_1^2 \eye$, we have  for $\forall t > 0$: 
\[
\alpha_{t} \leq (1 + 2\eta_t \lambda_1 + \eta_t^2 (\lambda_1^2 + V)) \alpha_{t-1}. 
\]
The result follows by using induction along with $\alpha_0=1$ and $1+ x \leq e^x$. 
\end{proof}

Using Lemma~\ref{lem:spectralboundE} we next bound the expected value of $\trace{\Vperp^\top \Bt{n} \BtTr{n} \Vperp}$. Ultimately this will allow us to bound the value $\trace{\Vperp^\top \Bt{n} \BtTr{n} \Vperp}$ with by Markov's inequality.

\begin{lemma}\label{lem:traceboundE}
For all $t \geq 0$ and $\eta_i \leq \frac{1}{\lambda_1}$ the following holds
\begin{multline*}
\E{\trace{\Vperp^\top \Bt{t} \BtTr{t} \Vperp}} \leq 
\exp\left(\sum_{j \in [t]} 2\eta_j \lambda_2 + \eta_j^2 \variancetotal \right) \cdot \left(d+\variance \sum_{i= 1}^{t}\eta_i^2 \exp
\left(\sum_{j\in[i]} 2\eta_j(\lambda_1-\lambda_2) \right)\right)
% \left[
% d \cdot \exp\left(3 \lambda_2 \sum_{j \in [t]} \eta_j \right)
% + \sum_{i \in [t]} \eta_i^2 \cdot \exp\left(
% 3\lambda_1 \sum_{j \in [i]} \eta_j + 3 \lambda_2 \sum_{j = i + 1}^{t} \eta_j
% \right)
% \right]
%\\
%& \enspace \enspace \enspace
%\cdot \exp\left(\sum_{i \in [t]} \eta_i^2 V\right)
~.
\end{multline*}
\end{lemma}
%
%\begin{proof}
%Let  $\alpha_t \eqdef \E{\trace{\UsperpTr \Bt{t} \BtTr{t} \Usperp}}$. We have
%\begin{align*}
%\alpha_t
%&= 
%\E{\trace{\Vperp^\top (\eye + \eta_t \A_t) \Bt{t - 1} \Bt{t - 1}^\top (\eye + \eta_t \A_t)^\top \Vperp}}
%\\
%&= 
%\E{\trace{
%		\Vperp^\top \Bt{t - 1} \Bt{t - 1}^\top \Vperp
%		+ 2\eta_t \Vperp^\top \Cov \Bt{t - 1} \Bt{t - 1}^\top \Vperp
%		+ \eta_t^2 \Vperp^\top \A_t \Bt{t - 1} \Bt{t - 1}^\top \A_t^\top \Vperp
%}}
%\\
%&\leq (1 + 2\eta_t \lambda_2) \alpha_{t - 1} + \eta_t^2 
%\E{\trace{\Vperp^\top \A_t \Bt{t - 1} \Bt{t - 1}^\top \A_t^\top \Vperp}}
%\end{align*}
%where in the last line we used the cyclic property of trace and definition of $\Vperp$ in
%\[
%\trace{\Vperp^\top \Cov \Bt{t - 1} \Bt{t - 1}^\top \Vperp}
%= \trace{\Bt{t - 1}^\top \Vperp \Vperp^\top \Cov \Bt{t - 1}}
%\leq \lambda_2 \trace{\Bt{t - 1}^\top \Vperp \Vperp^\top \Bt{t - 1}}
%= \alpha_{t - 1} ~.
%\]
%\end{proof}
\begin{proof}
Let $\alpha_t \eqdef \E{\trace{\Vperp^\top \Bt{t} \BtTr{t} \Vperp}}$. We first simplify $\alpha_t$ as follows:
\begin{align}
\alpha_t &= \iprod{\E{\Bt{t}\BtTr{t}}}{\Vperp \Vperp^\top }
= \iprod{\E{\Bt{t-1}\BtTr{t-1}}}{\E{\left(\eye+\eta_t\At{t}\right)
		\Vperp \Vperp^\top 
		\left(\eye+\eta_t\AtTrans{t}\right)} }.
\label{eqn:traceBoundE1}
\end{align}
Recall that $\E{\A_t}=\Cov$. Now, the second term on the right hand side can be bounded as follows:
\begin{align*}
&\E{\left(\eye+\eta_t\At{t}\right) \Vperp \Vperp^\top \left(\eye+\eta_t\AtTrans{t}\right)}, \\
&= \Vperp \Vperp^\top + \eta_t \Cov \Vperp \Vperp^\top + \eta_t \Vperp \Vperp^\top \Cov + \eta_t^2 \E{\At{t} \Vperp \Vperp^\top \AtTrans{t}}, \\
&= \Vperp \Vperp^\top + \eta_t \Cov \Vperp \Vperp^\top + \eta_t \Vperp \Vperp^\top \Cov + \eta_t^2 {\Cov \Vperp \Vperp^\top \Cov } + \eta_t^2 \E{\left(\At{t}-\Cov\right)\Vperp \Vperp^\top \trans{\left(\At{t}-\Cov\right)}}, \\
&\stackrel{\zeta_1}{\preceq} \Vperp \Vperp^\top + 2 \eta_t \lambda_2 \Vperp \Vperp^\top + \eta_t^2 \lambda_2^2 \Vperp \Vperp^\top+ \eta_t^2 \E{ \left(\At{t}-\Cov\right) \trans{\left(\At{t}-\Cov\right)} }, \\
&\stackrel{\zeta_2}{\preceq} \left(1+2 \eta_t \lambda_2 + \eta_t^2 \lambda_2^2 \right) \Vperp \Vperp^\top + \eta_t^2 \variance \eye = \left(1+2 \eta_t \lambda_2 + \eta_t^2 \lambda_2^2 + \eta_t^2 \variance \right) \Vperp \Vperp^\top + \eta_t^2 \variance \cdot \v_1 \v_1^\top,
\end{align*}
where $\zeta_1$ follows from the fact that $\Vperp$ is orthogonal to $\v_1$ and $\zeta_2$ follows from defintion of $\variance$. 

Plugging the above into~\eqref{eqn:traceBoundE1}, we get for all $t \geq 1$, 
\begin{align*}
\alpha_t &\leq \left(1+2 \eta_t \lambda_2 + \eta_t^2 (\lambda_2^2 + \variance)\right) \iprod{\E{\Bt{t-1}\BtTr{t-1}}}{\Vperp \Vperp^\top} + \eta_t^2 \variance \iprod{\E{\Bt{t-1}\BtTr{t-1}}}{\v_1 \v_1^\top}, \\
&\leq \left(1+2 \eta_t \lambda_2 + \eta_t^2 \variancetotal\right) \alpha_{t-1} + \eta_t^2 \variance \twonorm{\E{\Bt{t-1}\BtTr{t-1}}}, \\
&\leq \exp\left(2 \eta_t \lambda_{2} + \eta_t^2 \variancetotal \right) \alpha_{t-1} 
+ \eta_t^2 \variance \exp\left(
\sum_{i\in[t - 1]} \eta_i \lambda_1 + \eta_i^2 \variancetotal\right), %\mbox{(Using Lemma~\ref{lem:spectralboundE}).}
\end{align*}
where the last inequality follows from $1+x\leq e^x$ and using Lemma~\ref{lem:spectralboundE}. 

Recursing the above inequality, we obtain
\begin{align*}
\alpha_t 
&\leq
\sum_{i \in [t]}
%\left[
\eta_i^2 \variance
%\cdot 
\exp \left( 
\sum_{j = i + 1}^{t}
2\eta_j \lambda_2 + \eta_j^2 \variancetotal\right)
%\cdot  
\exp
\left(\sum_{j\in[i]} 2\eta_j \lambda_1 + \eta_j^2 \variancetotal \right)
%\right]
%\\
%& 
%\enspace \hspace*{3in} 
+ \exp\left(\sum_{j \in [t]} 2\eta_j \lambda_2 + \eta_j^2 \variancetotal\right) \alpha_0,\\
& \leq\exp\left(\sum_{j \in [t]} 2\eta_j \lambda_2 + \eta_j^2 \variancetotal\right) \left(\alpha_0+\variance \sum_{i =1}^{t}\eta_i^2 \exp
\left(\sum_{j\in[i]} 2\eta_j(\lambda_1-\lambda_2) + \eta_j^2 \variancetotal \right)\right)
\end{align*}
Since $\Bt{0} = \eye$ we see that $\alpha_0 = d - 1 \leq d$. Using that $\eta_i \leq \frac{1}{\lambda_1} \leq \frac{1}{\lambda_2}$ completes the proof.
\end{proof}

Next we provide the lemmas that will allow us to lower bound $\v_1^\top \Bt{t} \Bt{t}^\top \v_1$. In Lemma~\ref{lem:lowerboundE} we lower bound $\E{\v_1^\top \Bt{t} \Bt{t}^\top \v_1}$ and in Lemma~\ref{lem:lowerboundV} we upper bound $\Var{\v_1^\top \Bt{t} \Bt{t}^\top \v_1}$. Ultimately, the lower bound follows using Chebyshev's inequality. 

\begin{lemma}\label{lem:lowerboundE}
For all $t \geq 0$ and $\eta_i \geq 0$ we have
\[
\E{\v_1^\top \Bt{t} \Bt{t}^\top \v_1}
\geq \exp \left(\sum_{i \in [t]} 2\eta_i \lambda_1 - 4\eta_i^2 \lambda_1^2\right)
\]
If we further assume that $\eta_i \leq \frac{1}{4\cdot \max\{\lambda_1,M\}}$ then
$
\E{\v_1^\top \Bt{t} \Bt{t}^\top \v_1}
\geq \exp (\lambda_1 \sum_{i \in [t]} \eta_i)
$.
\end{lemma}
\begin{proof}
	Let $\beta_t \eqdef \E{\v_1^\top \Bt{t}\BtTr{t} \v_1}$. Since $\Bt{t} = \left(\eye+\eta_t \At{t}\right) \Bt{t-1}$, we can bound $\beta_t$ as
	\begin{align*}
		\beta_t 
		&= \iprod{\E{\Bt{t-1}\BtTr{t-1}}}{\E{\left(\eye+\eta_t \At{t}\right) \v_1 \v_1^\top \left(\eye+\eta_t \AtTrans{t}\right)}} \\
		&= \iprod{\E{\Bt{t-1}\BtTr{t-1}}}{\v_1 \v_1^\top + \eta_t \Cov \v_1 \v_1^\top + \eta_t \v_1 \v_1^\top \Cov + \eta_t^2 \E{\At{t} \v_1 \v_1^\top \ustarTr\AtTrans{t}}}
		\\
		&\geq \iprod{\E{\Bt{t-1}\BtTr{t-1}}}{\v_1 \v_1^\top + \lambda_1 \eta_t \v_1 \v_1^\top + \lambda_1 \eta_t \v_1 \v_1^\top} ~.
	\end{align*}
Consequently $\beta_t \geq (1 + 2\eta_t \lambda_1) \beta_{t - 1}$. Furthermore, $\Bt{0}=\eye$ and hence $\beta_0 = \twonorm{\v_1}^2 = 1$. Proceeding by induction and using that  $1 + x \geq \exp(x - x^2)$ for all $x \geq 0$ finishes the proof.
\end{proof}

\begin{lemma}\label{lem:lowerboundV}
For $t \geq 0$ suppose that $\eta_i \leq \frac{1}{4\cdot \max\{\lambda_1,M\}}$ for all $i \in [t]$ then.
\[
\E{\left(\v_1^\top \Bt{t}\BtTr{t} \v_1\right)^2}
\leq 
\exp\left( \sum_{i \in [t]} 4\eta_i \lambda_1 + 10\eta_i^2 \variancetotal \right)
\]
\end{lemma}
\begin{proof}
	Let $\Wts{t}{s} \eqdef \left(\eye + \eta_t \At{t}\right)\cdots \left(\eye + \eta_{t-s+1}\At{t-s+1}\right)$ and $\gamma_s \eqdef \E{\left(\v_1^\top \Wts{t}{s}\WtsTr{t}{s} \v_1\right)^2}$. Note that $\Wts{t}{t}=\Bt{t}$ and $\gamma_t=\E{\v_1^\top\Bt{t}\BtTr{t}\v_1}$. Now, 
	\begin{align}
		\gamma_t
%		&= \E{\ustarTr \Bt{t}\BtTr{t} \ustar \ustarTr \Bt{t}\BtTr{t} \ustar} \nonumber \\
%		&= \trace{\E{\BtTr{t} \ustar \ustarTr \Bt{t}\BtTr{t} \ustar \ustarTr \Bt{t}}} \nonumber \\
		&= \trace{\E{\WtsTr{t}{t}\v_1 \v_1^\top \Wts{t}{t} \WtsTr{t}{t}\v_1 \v_1^\top \Wts{t}{t}}} \nonumber \\
		&= \trace{\E{(\eye+\eta_1 \AtTrans{1}) \WtsTr{t}{t-1}\v_1 \v_1^\top \Wts{t}{t-1} (\eye+\eta_1 \At{1})(\eye+\eta_1 \AtTrans{1}) \WtsTr{t}{t-1} \v_1 \v_1^\top \Wts{t}{t-1} (\eye+\eta_1 \At{1})}} \nonumber \\
		&= \trace{\E{(\eye+\eta_1 \AtTrans{1}) \Gt{t-1} (\eye+\eta_1 \At{1}) (\eye+\eta_1 \AtTrans{1}) \Gt{t-1} (\eye+\eta_1 \At{1})}}, \label{eqn:upperboundV-1}
	\end{align}
	where  $ \Gt{t-1} \eqdef {\WtsTr{t}{t-1}\v_1 \v_1^\top \Wts{t}{t-1}}$. In order to bound the above quantity, we first bound the above expression for an arbitrary $\Gt{t-1} \equiv \G$. We then take an expectation over only $\At{1}$ and then finally take an expectation over $\Gt{t-1}$. That is, for an arbitrary fixed symmetric matrix $\G$, we have: 
	\begin{align}
		&\trace{\E{\left(\eye+\eta_1 \AtTrans{1}\right) \G \left(\eye+\eta_1 \At{1}\right)\left(\eye+\eta_1 \AtTrans{1}\right) \G \left(\eye+\eta_1 \At{1}\right)}} \nonumber \\
		&=\trace{\E{\left(\G + \eta_1 \At{1}^\top \G + \eta_1 \G \At{1} + \eta_1^2 \At{1}^\top \G \At{1}\right)^2}} \nonumber \\
		&= \mathrm{Tr}\left(\G^2 + \eta_1 \E{\AtTrans{1}} \G^2 + \eta_1 \G^2 \E{\At{1}} + \eta_1 \G\left(\E{ \At{1}} + \E{\AtTrans{1}}\right)\G \right. \nonumber \\
		& \quad + \eta_1^2 \E{\AtTrans{1} \G \At{1} \G} + \eta_1^2 \E{\AtTrans{1} \G \AtTrans{1} \G} + \eta_1^2 \E{\G \At{1} \G \At{1}} + \eta_1^2 \E{\G \AtTrans{1} \G \At{1}} \nonumber \\
		&\quad + \eta_1^2 \G \E{\At{1} \AtTrans{1}} \G 
		+ \eta_1^2 \E{\AtTrans{1} \G^2 \At{1}} + \eta_1^3 \E{\AtTrans{1} \G \left(\At{1}+\AtTrans{1}\right) \G \At{1}} 
		 \nonumber \\
		&\quad \left. + \eta_1^3 \E{\AtTrans{1} \G \At{1} \AtTrans{1} \G } + \eta_1^3 \E{ \G \At{1}\AtTrans{1} \G \At{1}} + \eta_1^4 \E{\AtTrans{1} \G \At{1}\AtTrans{1} \G \At{1}}   \right) \nonumber \\
		&= \trace{\G^2} + 4 \eta_1 \trace{\Cov \G^2} + 2 \eta_1^2 \trace{\E{\At{1} \AtTrans{1}} \G^2} + \eta_1^2 \trace{\E{\AtTrans{1} \G \At{1} \G}} \nonumber \\
		& \quad + \eta_1^2 \trace{\E{\AtTrans{1} \G \AtTrans{1} \G}} + \eta_1^2 \trace{\E{\G \At{1} \G \At{1}}} + \eta_1^2 \trace{\E{\G \AtTrans{1} \G \At{1}}} \nonumber \\
		&\quad + 2 \eta_1^3 \trace{\E{\AtTrans{1} \G \left(\At{1}+\AtTrans{1}\right) \G \At{1}}} + \eta_1^4 \trace{\E{\AtTrans{1} \G \At{1}\AtTrans{1} \G \At{1}}} \label{eqn:upperboundV-2}
%		&\leq \trace{\G^2} + 4 \eta_t \trace{\Cov \G^2} + 4 \eta_t^2 \trace{\E{\At{1} \G^2 \At{1}}} + \eta_t^2 \trace{\E{\At{1}^2} \G^2} \nonumber \\
%		&\quad  + 4 \eta_t^3 \trace{\E{\At{1} \G \At{1}^2 \G }} + \eta_t^4 \trace{\E{\At{1} \G \At{1}^2 \G \At{1}}} \nonumber \\
%		&= \trace{\G^2} + 4 \eta_t \trace{\Cov \G^2} + 5 \eta_t^2 \trace{\E{\At{1}^2} \G^2} 
%		+ 4 \eta_t^3 \trace{\E{\At{1} \G \At{1} \G \At{1}}} \nonumber \\
%		&+ \eta_t^4 \trace{\E{\At{1} \G \At{1}^2 \G \At{1}}},
	\end{align}
	We now bound the various terms above as follows. Each of the second order terms can be bounded using Lemma~\ref{lem:trace-cauchy-schwarz} as follows:
	\begin{align}
		\E{\trace{\AtTrans{1} \G \At{1} \G}} &\leq \frac{1}{2} \E{\frob{\AtTrans{1} \G}^2 + \frob{\At{1}\G}^2} \nonumber \\
		&= \frac{1}{2} \left(\trace{\G\E{\At{1}\AtTrans{1}}\G + \G\E{\AtTrans{1}\At{1}}\G} \right) \leq  (\variance + \lambda_1^2 ) \trace{\G^2}. \label{eqn:secondorder}
	\end{align} 
	The third order terms can be bounded as follows:
	\begin{align}
		\E{\trace{\AtTrans{1}\G\At{1}\G\At{1}}} &\leq \E{\twonorm{\At{1}} \trace{\AtTrans{1}\G \G\At{1}}} \nonumber \\
		&\leq \left(\upperbound + \lambda_1\right) \trace{\G\E{\At{1}\AtTrans{1}}\G}
		\leq \left(\upperbound + \lambda_1\right) \variancetotal \cdot
		 \trace{\G^2}. \label{eqn:thirdorder}
	\end{align}
	where we used the assumption that $\twonorm{\At{1}} \leq \twonorm{\At{1}-\Cov} + \twonorm{\Cov} \leq \upperbound + \lambda_1$ with probability $1$. Finally the fourth order term can be bounded as
	\begin{align}
		\trace{\E{\AtTrans{1} \G \At{1} \AtTrans{1} \G \At{1}}}
		&\leq \left(\upperbound + \lambda_1\right)^2 \trace{\G^2 \E{\At{1} \AtTrans{1}}}
		\leq \left(\upperbound + \lambda_1\right)^2 \variancetotal \cdot \trace{\G^2}. \label{eqn:fourthorder}
	\end{align}
	Plugging~\eqref{eqn:secondorder},~\eqref{eqn:thirdorder} and~\eqref{eqn:fourthorder} into~\eqref{eqn:upperboundV-2} tells us that
	\begin{align*}
		&\trace{\E{\left(\eye+\eta_1 \AtTrans{1}\right) \G \left(\eye+\eta_1 \At{1}\right)\left(\eye+\eta_1 \AtTrans{1}\right) \G \left(\eye+\eta_t \At{1}\right)}} \\
		&\quad \leq \trace{\G^2} + 4 \eta_1 \lambda_1 \trace{ \G^2} + 5 \eta_1^2 \variancetotal \cdot \trace{\G^2}  \\
		& \qquad + 4 \eta_1^3 \left(\upperbound + \lambda_1 \right) \variancetotal \cdot \trace{\G^2} + \eta_1^4 \left(\upperbound + \lambda_1 \right)^2 \variancetotal \cdot \trace{\G^2} \\
		&\quad = \left( 1 + 4\eta_1 \lambda_1 + 5\eta_1^2 \variancetotal + 4 \eta_1^3 \left(\upperbound + \lambda_1\right) \variancetotal + \eta_1^4 \left(\upperbound + \lambda_1\right)^2 \variancetotal \right) \trace{\G^2} \\
		&\quad \leq \exp\left( 4\eta_1 \lambda_1 + 10\eta_1^2 \variancetotal \right) \trace{\G^2}
	\end{align*}
	where in the last line we used that $\eta_i \leq \frac{1}{4 \max\{\upperbound,\lambda_1\}}$ and that $1 + x \leq \exp(x)$

	Using the value $\G = \Gt{t-1}={\WtsTr{t}{t-1}\v_1 \v_1^\top \Wts{t}{t-1}}$ and plugging the above into~\eqref{eqn:upperboundV-1}, we have
	\begin{align*}
		\gamma_t &= \trace{\E{\left(\eye+\eta_1 \AtTrans{1}\right) \Gt{t-1} \left(\eye+\eta_1 \At{1}\right) \left(\eye+\eta_1 \AtTrans{1}\right) \Gt{t-1} \left(\eye+\eta_1 \At{1}\right)}} \\
		&\leq \exp\left( 4\eta_1 \lambda_1 + 10\eta_1^2 \variancetotal  \right) \E{\trace{\Gt{t-1}^2}} 
		%\\
		%&
		= \exp\left( 4\eta_1 \lambda_1 + 10\eta_1^2 \variancetotal \right) \gamma_{t-1},
	\end{align*}
	where we used the fact that $\gamma_{t-1} = \E{\trace{\Gt{t-1}^2}}$. Since $\gamma_0 = 1$, induction proves the lemma.
\end{proof}


We now have everything to prove Theorem~\ref{thm:main1}.

\begin{proof}[Proof of Theorem~\ref{thm:main1}]
	As discussed in Section~\ref{sec:approach} the main idea of this proof to use that Algorithm~\ref{algo:oja} is essentially one step of power method for the matrix $\Bt{n}$ and use Lemma~\ref{lem:onesteppowermethod} to bound the error.  To this end, we lower and upper bound $\v_1^\top\Bt{n}\BtTr{n}\v_1$ and	$\trace{\Vperp^\top\Bt{n}\BtTr{n}\Vperp}$, respectively.

First, using Chebyshev's inequality, we have: 
	\begin{align*}
		\P{\abs{\v_1^\top\Bt{n}\BtTr{n}\v_1 - \E{\v_1^\top\Bt{n}\BtTr{n}\v_1}} > \frac{1}{\sqrt{\delta}}\sqrt{\Var{\v_1^\top\Bt{n}\BtTr{n}\v_1}}} < \delta.
	\end{align*}
	So with probability greater than $1-\delta$, the following holds: 
	\begin{align}
		&\v_1^\top\Bt{n}\BtTr{n}\v_1 > \E{\v_1^\top\Bt{n}\BtTr{n}\v_1} - \frac{1}{\sqrt{\delta}}\sqrt{\Var{\v_1^\top\Bt{n}\BtTr{n}\v_1}} \notag\\
		&= \E{\v_1^\top\Bt{n}\BtTr{n}\v_1} - \frac{1}{\sqrt{\delta}}\sqrt{\E{\left(\v_1^\top\Bt{n}\BtTr{n}\v_1\right)^2} - \E{\v_1^\top\Bt{n}\BtTr{n}\v_1}^2}  \notag\\
		&\stackrel{\zeta_1}{\geq} \exp\left(2\lambda_1 \sum_{i=1}^{n}\eta_i -4\lambda_1^2\sum_{i=1}^{n}\eta_i^2\right) \times \left( 1 - \frac{1}{\sqrt{\delta}}\sqrt{ \exp\left(18\sum_{i=1}^{n} \eta_i^2  
			\variancetotal \right) -  1}\right) \notag \\
\label{eqn:main1}
		\end{align}
where $\zeta_1$ follows from Lemma~\ref{lem:lowerboundE} and \ref{lem:lowerboundV}. 

		Furthermore, using Lemma~\ref{lem:traceboundE} and Markov's inequality, we have with probability at least $ 1-\delta$,
		\begin{multline}
\trace{\Vperp^\top \Bt{t} \BtTr{t} \Vperp}
\leq \frac{\exp\left(\sum_{i \in [n]} 2\eta_i \lambda_2 + \eta_i^2 \variancetotal \right)}{\delta} \cdot \left(d+\variance \sum_{i=1}^{n}\eta_i^2 \exp
\left(\sum_{j\in[i]} 2\eta_j(\lambda_1-\lambda_2) \right)\right)
% \frac{\exp(\sum_{i \in [t]} \eta_i^2 \variance)}{\delta}\cdot 
% \left[
% d \cdot \exp(3 \lambda_2 \sum_{j \in [t]} \eta_j )
% + \sum_{i \in [t]} \eta_i^2 \cdot \exp(
% 3\lambda_1 \sum_{j \in [i]} \eta_j + 3 \lambda_2 \sum_{j = i + 1}^{t} \eta_j)
% \right]%\\ \cdot \exp\left(\sum_{i \in [t]} \eta_i^2 V\right)
~. \label{eqn:main2}
		\end{multline}
		Consequently with probability at least $1 - 2\delta$ both  \eqref{eqn:main1} and \eqref{eqn:main2} hold and therefore the result follows by  Lemma~\ref{lem:onesteppowermethod} and choosing a $\delta$ that is smaller by a constant.
%		\begin{align}
%			&1-(\wt{n}^\top \v_1)^2 \leq \frac{\trace{\Vperp^\top \Bt{n} \BtTr{n} \Vperp}}{\v_1^\top\Bt{n}\BtTr{n}\v_1}\notag \\&\leq \frac{1}{\delta\cdot Q} \cdot \left(d\cdot \exp\left(-2(\lambda_1-\lambda_2)\sum_{i \in [n]} \eta_i  + 5 \overline{\variance}\sum_{i\in[n]}\eta_i^2\right)+\variance \sum_{i=1}^{n}\eta_i^2 \exp
%(-\sum_{j=i+1}^n 2\eta_j(\lambda_1-\lambda_2) )\right),\label{eqn:main3}
%		\end{align}
%where $Q= \left( 1 - \frac{1}{\sqrt{\delta}}\sqrt{ \exp\left(18\overline{\variance}\sum_{i=1}^{n} \eta_i^2 \right) -  1}\right)$.
\end{proof}

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

notations.tex
%% General machinery
\newcommand{\be}{\begin{equation}}
\newcommand{\ee}{\end{equation}}
\newcommand{\argmin}{\mathop{\rm argmin}}
\newcommand{\argmax}{\mathop{\rm argmax}}
\newcommand{\AutoAdjust}[3]{\mathchoice{ \left #1 #2  \right #3}{#1 #2 #3}{#1 #2 #3}{#1 #2 #3} }
\newcommand{\super}[1]{^{(#1)}}
\newcommand{\dd}{\mathrm{d}}
\newcommand{\vect}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\mat}[1]{\ensuremath{\mathbf{#1}}}
\newcommand{\grad}{\bigtriangledown}
\newcommand{\hess}{\bigtriangledown^2}
%%%%%%%%%%%%% New macros %%%%%%%%%%%%%%%%%%
\newcommand{\twonorm}[1]{\left\| {#1} \right\|_2}
\newcommand{\trnorm}[1]{\left\| {#1} \right\|_{\ast}}
\newcommand{\iprod}[2]{\left\langle #1, #2 \right\rangle}
\newcommand{\frob}[1]{\left\|#1\right\|_F}
%\newcommand{\set}[1]{\left\{#1\right\}}
\newcommand{\order}[1]{\ensuremath{\mathcal{O}\left(#1\right)}}
\newcommand{\otilde}[1]{\ensuremath{\widetilde{\mathcal{O}}\left(#1\right)}}
\newcommand{\Om}[1]{\ensuremath{\Omega\left(#1\right)}}
\newcommand{\distop}[2]{\ensuremath{dist\left(#1,#2\right)}}
\newcommand{\eqdef}{\stackrel{\triangle}{=}}
\newcommand{\defeq}{\stackrel{\triangle}{=}}
\newcommand{\abs}[1]{\left|#1\right|}
\newcommand{\infnorm}[1]{\left\|#1\right\|_{\infty}}
\newcommand{\trace}[1]{\mathrm{Tr}\left(#1\right)}
\renewcommand{\det}[1]{{\rm det}\left(#1\right)}
\newcommand{\logdet}[1]{{\rm logdet}\left(#1\right)}
\newcommand{\inv}[1]{{#1}^{-1}}
\newcommand{\trans}[1]{{#1}^{\top}}
\newcommand{\wh}[1]{\widehat{#1}}
%\newcommand{\citep}{\cite}
%\newcommand{\citet}{\cite}
%\newcommand{\citealp}{\cite}
\def\nn{\nonumber}

%%%%%from other papers
\newcommand{\E}[1]{\mathbb{E}\left[#1\right]}
\newcommand{\Eop}{\mathbb{E}}
\renewcommand{\P}[1]{\mathbb{P}\left[#1\right]}
\newcommand{\Var}[1]{\mathrm{Var}\left[#1\right]}
\newcommand\R{\mathbb{R}}
\newcommand\N{\mathbb{N}}
\renewcommand\S{\mathbb{S}}
\def\Pbb{\mathbb{P}}
\newcommand\ch{\widehat{c}}
\newcommand\ct{\widetilde{c}}

\newcommand\Rn{\mathbb{R}^n}
\newcommand\Rnn{\mathbb{R}^{n\times n}}

\newcommand{\singmin}[1]{\ensuremath{\sigma_{\min}\left(#1\right)}}
\newcommand{\singmax}[1]{\ensuremath{\sigma_{\max}\left(#1\right)}}
\newcommand{\singi}[1]{\ensuremath{\sigma_{i}\left(#1\right)}}
\newcommand{\singj}[1]{\ensuremath{\sigma_{#1}}}
\newcommand{\sigstar}[1]{\ensuremath{\sigma_{#1}^*}}
\newcommand{\lammin}[1]{\ensuremath{\lambda_{\min}\left(#1\right)}}
\newcommand{\lami}[1]{\ensuremath{\lambda_{i}\left(#1\right)}}

\newcommand\eps{\varepsilon}

\newcommand{\ind}[1]{\mathds{1}\left[\vphantom{\sum}#1\right]}
\newcommand{\phiu}[1]{\phi_{\u}\left(#1\right)}
\newcommand{\phiustar}[1]{\phi_{\ustar}\left(#1\right)}
\newcommand{\phiuone}[1]{\phi_{\u_1}\left(#1\right)}
\newcommand{\phiv}[1]{\phi_{\v}\left(#1\right)}
\newcommand{\poly}[1]{\textrm{poly}\left(#1\right)}
\newcommand{\polylog}[1]{\textrm{polylog}\left(#1\right)}

\newcommand{\M}{\mat{M}}
\newcommand\Mmh{\mat{M}^{-1/2}}
\newcommand{\Cov}{\mat{\Sigma}}
\newcommand{\A}{\mat{A}}
\newcommand{\At}[1]{\mat{A}_{#1}}
\newcommand{\AtTrans}[1]{\mat{A}_{{#1}}^{\top}}
\newcommand{\B}{\mat{B}}
\newcommand{\Bt}[1]{\mat{B}_{#1}}
\newcommand{\BtTr}[1]{\mathbf{B}_{#1}^{\top}}
\newcommand{\C}{\mat{C}}
\newcommand{\D}{\mat{D}}
\newcommand{\G}{\mat{G}}
\newcommand{\W}{\mat{W}}
\newcommand{\V}{\mat{V}}
\newcommand{\Z}{\mat{Z}}
\newcommand{\mSigma}{\mat{\Sigma}}
\newcommand{\mZero}{\mat{0}}
\newcommand{\Wts}[2]{\mat{W}_{#1,#2}}
\newcommand{\Gt}[1]{{\mat{G}_{#1}}}
\newcommand{\WtsTr}[2]{\mathbf{W}_{#1,#2}^{\top}}
\newcommand{\Et}[1][t]{\mat{E_{#1}}}
\newcommand{\Etp}{\Et[t+1]}
\newcommand{\Errt}[1][t]{\mat{\bigtriangleup_{#1}}}
\newcommand\cnM{\kappa}
\newcommand{\cn}[1]{\kappa\left(#1\right)}
\newcommand\X{\mat{X}}
\newcommand\fstar{f_*}
\newcommand\Xt[1]{\mat{X_{#1}}}
\newcommand\ut[1][t]{{u_{#1}}}
\newcommand\Xtinv{\inv{\Xt}}
\newcommand\Xtp{\mat{X_{t+1}}}
\newcommand\Xtpinv{\inv{\left(\mat{X_{t+1}}\right)}}
\newcommand\U{\mat{U}}
\newcommand\UTr{\trans{\mat{U}}}
\newcommand{\Ut}[1][t]{\mat{U_{#1}}}
\newcommand{\Utinv}{\inv{\Ut}}
\newcommand{\UtTr}[1][t]{\trans{\mat{U_{#1}}}}
\newcommand\Utp{\mat{U_{t+1}}}
\newcommand\UtpTr{\trans{\mat{U}_{t+1}}}
\newcommand\Utptild{\mat{\widetilde{U}_{t+1}}}
\newcommand\Us{\mat{U^*}}
\newcommand{\Usperp}{\mat{U}^*_{\perp}}
\newcommand{\UsperpTr}{\trans{\Usperp}}
\newcommand\UsTr{\trans{\mat{U^*}}}
\newcommand{\Vperp}{\mat{V}_{\perp}}
\newcommand{\Sigs}{\mat{\Sigma}}
\newcommand{\Sigsmh}{\Sigs^{-1/2}}
\newcommand{\eye}{\mat{I}}
\newcommand{\twonormbound}{\left(4+\DPhi{\M}{\Xt[0]}\right)\twonorm{\M}}
\newcommand{\lamj}{\lambda_j}

\newcommand{\vbar}{\overline{v}}
\newcommand{\Mbar}{\overline{M}}
\renewcommand\u{\vect{u}}
\newcommand\x{\vect{x}}
\newcommand\e{\vect{e}}
\newcommand\uTr{\trans{\vect{u}}}
\newcommand\ustar{{\vect{u}^*}}
\newcommand\ustarTr{{{\vect{u}^*}^{\top}}}
\newcommand\utilde{{\widetilde{\vect{u}}}}
\newcommand\vtilde{{\widetilde{\vect{v}}}}
\newcommand{\valpha}{\vect{\alpha}}
\newcommand\vg{\vect{g}}
\newcommand\Utildeperp{{\widetilde{\mat{U}}}_{\perp}}
\newcommand\UtildeperpTr{{\widetilde{\mat{U}}}_{\perp}^{\top}}
\newcommand\Vtildeperp{{\widetilde{\mat{V}}}_{\perp}}
\newcommand\VtildeperpTr{{\widetilde{\mat{V}}}_{\perp}^{\top}}
\renewcommand\v{\vect{v}}
\newcommand{\va}{\vect{a}}
\newcommand\vTr{\trans{\vect{v}}}
\newcommand\w{\vect{w}}
\newcommand{\wt}[1]{\vect{w_{#1}}}
\newcommand\wTr{\trans{\vect{w}}}
\newcommand\wperp{\vect{w}_{\perp}}
\newcommand\wperpTr{\trans{\vect{w}_{\perp}}}
\newcommand\wj{\vect{w_j}}
\newcommand\vj{\vect{v_j}}
\newcommand\wjTr{\trans{\vect{w_j}}}
\newcommand\vjTr{\trans{\vect{v_j}}}

\newcommand{\DPhi}[2]{\ensuremath{D_{\Phi}\left(#1,#2\right)}}
\newcommand\matmult{{\omega}}
\newcommand\tzero{{t_0}}

% Objects
\newcommand{\dist}{\mathcal{D}}

% Our parameters
\newcommand{\variance}{\mathcal{V}}
\newcommand{\variancetotal}{\overline{\mathcal{V}}}
\newcommand{\upperbound}{\mathcal{M}}
\newcommand{\ninit}{n_{\textrm{init}}}
\newcommand{\nfinal}{n_{\textrm{final}}}

Ojasalgo.tex
\documentclass[11pt]{article}

\usepackage{amsmath,amsbsy,amsfonts,amssymb,amsthm,fullpage,units,url}
% The following packages will be automatically loaded:
% amsmath, amssymb, natbib, graphicx, url, algorithm2e
\usepackage{times,caption}


%\documentclass[11pt]{article}

\usepackage{units,stackrel}
%\usepackage{graphicx}
\usepackage{algorithmic,algorithm}
\usepackage{color,cases}
\usepackage{hyperref}
\usepackage{dsfont}
\usepackage[numbers]{natbib}
%\usepackage[lined,boxed,commentsnumbered]{algorithm2e}
%\usepackage{array,pifont}
\usepackage{array}

\newcolumntype{C}[1]{>{\centering\let\newline\\\arraybackslash\hspace{0pt}}m{#1}}
\newtheorem{claim}{Claim}
\newtheorem{corollary}{Corollary}[section]
\newtheorem{assumption}{Assumption}[section]
\newtheorem{lemma}{Lemma}[section]
\newtheorem{proposition}{Proposition}[section]
\newtheorem{theorem}{Theorem}[section]
\newtheorem{fact}{Fact}[section]
\theoremstyle{definition}
\newtheorem{condition}{Condition}[section]
\newtheorem{definition}{Definition}
\newtheorem{example}{Example}
%\theoremstyle{remark}
\newtheorem{remark}{Remark}

\newcommand{\sidford}{\textcolor[rgb]{0,0.5,0}{\textbf{Sidford}: }\textcolor[rgb]{0,0.5,0}}
\newcommand{\Chi}{\textcolor[rgb]{0,0.5,0.5}{\textbf{Chi}: }\textcolor[rgb]{0,0.5,0.5}}

\include{notations}


\newcommand{\todo}[1]{{\color{red}{TODO:#1}}}
\definecolor{WildStrawberry}{RGB}{255,67,164}

\newcommand{\praneeth}[1]{{\color{blue}{PN: #1}}}
\newcommand{\pj}[1]{{\color{magenta}{Prateek: #1}}}
\renewcommand{\cite}{~\citep}


\title{Streaming PCA: Matching Matrix Bernstein and Near-Optimal Finite Sample Guarantees for Oja's Algorithm}
%\date{}
%\Runningtitle{Near-Optimal Finite Sample Guarantees for Oja's Algorithm}
%\ShortHeadings{short title}{short authors}
\author{Prateek Jain\footnote{Microsoft Research India. Email: prajain@microsoft.com} \and Chi Jin\footnote{UC Berkeley. Email: chijin@cs.berkeley.edu} \and Sham M. Kakade\footnote{University of Washington. Email: sham@cs.washington.edu} \and Praneeth Netrapalli\footnote{Microsoft Research New England. Email: praneeth@microsoft.com} \and Aaron Sidford\footnote{Microsoft Research New England. Email: asid@microsoft.com}}


\begin{document}
\maketitle
%\ShortHeadings{short title}{short authors}

\begin{abstract}
This work provides improved guarantees for
% solving a key problem in data analysis, 
streaming principle component analysis (PCA). Given
$\A_1, \ldots, \A_n\in \R^{d\times d}$ sampled independently from
distributions satisfying $\E{\A_i} = \Cov$ for $\Cov \succeq \mZero$,
this work provides an $O(d)$-space linear-time single-pass streaming algorithm
for estimating the top eigenvector of $\Cov$. The algorithm nearly
matches (and in certain cases improves upon) the accuracy obtained by
the standard batch method that computes top eigenvector of the
empirical covariance $\frac{1}{n} \sum_{i \in [n]} \A_i$ as analyzed
by the matrix Bernstein inequality. Moreover, to achieve constant
accuracy, our algorithm improves upon the best previous known sample
complexities of streaming algorithms by either a multiplicative factor
of $O(d)$ or $1/\mathrm{gap}$ where $\mathrm{gap}$ is the relative
distance between the top two eigenvalues of $\Cov$.


These results are achieved through a novel analysis of the classic Oja's
algorithm, one of the oldest and most popular algorithms for
streaming PCA. In particular, this work shows that simply picking a random initial point
$\w_0$ and applying the update rule
$\w_{i + 1} = \w_i + \eta_i \A_i \w_i$ suffices to accurately estimate
the top eigenvector, with a suitable choice of
$\eta_i$. We believe our result sheds light on how to efficiently
perform streaming PCA both in theory and in practice, and we hope that
our analysis may serve as the basis for analyzing many variants and
extensions of streaming PCA.
\end{abstract}
\newpage
\input{intro}
\input{related}

\input{prelims}

\input{approach}

\input{results}
\input{lemmas}

\input{conclusion}

\section{Acknowledgements}
Sham Kakade acknowledges funding from the Washington Research Foundation for innovation in
Data-intensive Discovery.

\newpage

%\bibliography{bib}
%\bibliographystyle{IEEEtranN}
% Generated by IEEEtranN.bst, version: 1.13 (2008/09/30)
\begin{thebibliography}{27}
	\providecommand{\natexlab}[1]{#1}
	\providecommand{\url}[1]{#1}
	\csname url@samestyle\endcsname
	\providecommand{\newblock}{\relax}
	\providecommand{\bibinfo}[2]{#2}
	\providecommand{\BIBentrySTDinterwordspacing}{\spaceskip=0pt\relax}
	\providecommand{\BIBentryALTinterwordstretchfactor}{4}
	\providecommand{\BIBentryALTinterwordspacing}{\spaceskip=\fontdimen2\font plus
		\BIBentryALTinterwordstretchfactor\fontdimen3\font minus
		\fontdimen4\font\relax}
	\providecommand{\BIBforeignlanguage}[2]{{%
			\expandafter\ifx\csname l@#1\endcsname\relax
			\typeout{** WARNING: IEEEtranN.bst: No hyphenation pattern has been}%
			\typeout{** loaded for the language `#1'. Using the pattern for}%
			\typeout{** the default language instead.}%
			\else
			\language=\csname l@#1\endcsname
			\fi
			#2}}
	\providecommand{\BIBdecl}{\relax}
	\BIBdecl
	
	\bibitem[Jolliffe(2002)]{jolliffe2002principal}
	I.~Jolliffe, \emph{Principal component analysis}.\hskip 1em plus 0.5em minus
	0.4em\relax Wiley Online Library, 2002.
	
	\bibitem[Hall et~al.(1998)Hall, Marshall, and Martin]{hall1998incremental}
	P.~M. Hall, A.~D. Marshall, and R.~R. Martin, ``Incremental eigenanalysis for
	classification.'' in \emph{BMVC}, vol.~98.\hskip 1em plus 0.5em minus
	0.4em\relax Citeseer, 1998, pp. 286--295.
	
	\bibitem[Weng et~al.(2003)Weng, Zhang, and Hwang]{weng2003candid}
	J.~Weng, Y.~Zhang, and W.-S. Hwang, ``Candid covariance-free incremental
	principal component analysis,'' \emph{Pattern Analysis and Machine
		Intelligence, IEEE Transactions on}, vol.~25, no.~8, pp. 1034--1040, 2003.
	
	\bibitem[Ross et~al.(2008)Ross, Lim, Lin, and Yang]{ross2008incremental}
	D.~A. Ross, J.~Lim, R.-S. Lin, and M.-H. Yang, ``Incremental learning for
	robust visual tracking,'' \emph{International Journal of Computer Vision},
	vol.~77, no. 1-3, pp. 125--141, 2008.
	
	\bibitem[Krasulina(1970)]{krasulina1970method}
	T.~Krasulina, ``Method of stochastic approximation in the determination of the
	largest eigenvalue of the mathematical expectation of random matrices,''
	\emph{Automatation and Remote Control}, pp. 50--56, 1970.
	
	\bibitem[Oja(1982)]{oja1982simplified}
	E.~Oja, ``Simplified neuron model as a principal component analyzer,''
	\emph{Journal of mathematical biology}, vol.~15, no.~3, pp. 267--273, 1982.
	
	\bibitem[Balsubramani et~al.(2013)Balsubramani, Dasgupta, and
	Freund]{BalsubramaniDF13}
	A.~Balsubramani, S.~Dasgupta, and Y.~Freund, ``The fast convergence of
	incremental pca,'' in \emph{Advances in Neural Information Processing
		Systems}, 2013, pp. 3174--3182.
	
	\bibitem[Mitliagkas et~al.(2013)Mitliagkas, Caramanis, and
	Jain]{MitgliakisCJ13}
	I.~Mitliagkas, C.~Caramanis, and P.~Jain, ``Memory limited, streaming pca,'' in
	\emph{Advances in Neural Information Processing Systems}, 2013, pp.
	2886--2894.
	
	\bibitem[Sa et~al.(2015)Sa, Re, and Olukotun]{SaOR15}
	C.~D. Sa, C.~Re, and K.~Olukotun, ``Global convergence of stochastic gradient
	descent for some non-convex matrix problems,'' in \emph{Proceedings of the
		32nd International Conference on Machine Learning, {ICML} 2015, Lille,
		France, 6-11 July 2015}, 2015, pp. 2332--2341.
	
	\bibitem[Vershynin(2010)]{vershynin2010introduction}
	R.~Vershynin, ``Introduction to the non-asymptotic analysis of random
	matrices,'' \emph{arXiv preprint arXiv:1011.3027}, 2010.
	
	\bibitem[Tropp(2012)]{Tropp12}
	J.~A. Tropp, ``User-friendly tail bounds for sums of random matrices,''
	\emph{Foundations of Computational Mathematics}, vol.~12, no.~4, pp.
	389--434, 2012.
	
	\bibitem[Johnstone(2001)]{johnstone2001distribution}
	I.~M. Johnstone, ``On the distribution of the largest eigenvalue in principal
	components analysis,'' \emph{Annals of statistics}, pp. 295--327, 2001.
	
	\bibitem[Wedin(1972)]{wedin1972perturbation}
	P.-{\AA}. Wedin, ``Perturbation bounds in connection with singular value
	decomposition,'' \emph{BIT Numerical Mathematics}, vol.~12, no.~1, pp.
	99--111, 1972.
	
	\bibitem[Jin et~al.(2015)Jin, Kakade, Musco, Netrapalli, and
	Sidford]{jin2015robust}
	C.~Jin, S.~M. Kakade, C.~Musco, P.~Netrapalli, and A.~Sidford, ``Robust
	shift-and-invert preconditioning: Faster and more sample efficient algorithms
	for eigenvector computation,'' \emph{arXiv preprint arXiv:1510.08896}, 2015.
	
	\bibitem[Golub and Van~Loan(2012)]{golub2012matrix}
	G.~H. Golub and C.~F. Van~Loan, \emph{Matrix computations}.\hskip 1em plus
	0.5em minus 0.4em\relax JHU Press, 2012, vol.~3.
	
	\bibitem[Hardt and Price(2014)]{HardtP14}
	M.~Hardt and E.~Price, ``The noisy power method: {A} meta algorithm with
	applications,'' in \emph{Advances in Neural Information Processing Systems
		27: Annual Conference on Neural Information Processing Systems 2014, December
		8-13 2014, Montreal, Quebec, Canada}, 2014, pp. 2861--2869.
	
	\bibitem[Cohen et~al.(2010)Cohen, Lee, Miller, Pachocki, and
	Sidford]{cohen2016median}
	M.~Cohen, Y.~T. Lee, G.~Miller, J.~Pachocki, and A.~Sidford, ``Geometric median
	in nearly linear time,'' \emph{To Appear in 48th Annual Symposium on the
		Theory of Computing (STOC) 2016}, 2010.
	
	\bibitem[Clarkson and Woodruff(2009)]{ClarksonW09}
	K.~L. Clarkson and D.~P. Woodruff, ``Numerical linear algebra in the streaming
	model,'' in \emph{Proceedings of the forty-first annual ACM symposium on
		Theory of computing}.\hskip 1em plus 0.5em minus 0.4em\relax ACM, 2009, pp.
	205--214.
	
	\bibitem[Liberty(2013)]{Liberty13}
	E.~Liberty, ``Simple and deterministic matrix sketching,'' in \emph{Proceedings
		of the 19th ACM SIGKDD international conference on Knowledge discovery and
		data mining}.\hskip 1em plus 0.5em minus 0.4em\relax ACM, 2013, pp. 581--588.
	
	\bibitem[Nelson and Nguy{\^e}n(2013)]{nelson2013osnap}
	J.~Nelson and H.~L. Nguy{\^e}n, ``Osnap: Faster numerical linear algebra
	algorithms via sparser subspace embeddings,'' in \emph{Foundations of
		Computer Science (FOCS), 2013 IEEE 54th Annual Symposium on}.\hskip 1em plus
	0.5em minus 0.4em\relax IEEE, 2013, pp. 117--126.
	
	\bibitem[Cohen et~al.(2015)Cohen, Nelson, and Woodruff]{cohen2015optimal}
	M.~B. Cohen, J.~Nelson, and D.~P. Woodruff, ``Optimal approximate matrix
	product in terms of stable rank,'' \emph{arXiv preprint arXiv:1507.02268},
	2015.
	
	\bibitem[Ghashami et~al.(2015)Ghashami, Liberty, Phillips, and
	Woodruff]{GhashamiLPW15}
	M.~Ghashami, E.~Liberty, J.~M. Phillips, and D.~P. Woodruff, ``Frequent
	directions: Simple and deterministic matrix sketching,'' \emph{arXiv preprint
		arXiv:1501.01711}, 2015.
	
	\bibitem[Boutsidis et~al.(2015)Boutsidis, Garber, Karnin, and
	Liberty]{BoutsidisGKL15}
	C.~Boutsidis, D.~Garber, Z.~Karnin, and E.~Liberty, ``Online principal
	components analysis,'' in \emph{Proceedings of the Twenty-Sixth Annual
		ACM-SIAM Symposium on Discrete Algorithms}.\hskip 1em plus 0.5em minus
	0.4em\relax SIAM, 2015, pp. 887--901.
	
	\bibitem[Shamir(2015)]{Shamir15}
	O.~Shamir, ``A stochastic {PCA} and {SVD} algorithm with an exponential
	convergence rate,'' in \emph{Proceedings of the 32nd International Conference
		on Machine Learning, {ICML} 2015, Lille, France, 6-11 July 2015}, 2015, pp.
	144--152.
	
	\bibitem[Garber and Hazan(2015)]{garber2015fast}
	D.~Garber and E.~Hazan, ``Fast and simple pca via convex optimization,''
	\emph{arXiv preprint arXiv:1509.05647}, 2015.
	
	\bibitem[Warmuth and Kuzmin(2006)]{WarmuthK06}
	M.~K. Warmuth and D.~Kuzmin, ``Randomized {PCA} algorithms with regret bounds
	that are logarithmic in the dimension,'' in \emph{Advances in Neural
		Information Processing Systems 19, Proceedings of the Twentieth Annual
		Conference on Neural Information Processing Systems, Vancouver, British
		Columbia, Canada, December 4-7, 2006}, 2006, pp. 1481--1488.
	
	\bibitem[Garber et~al.(2015)Garber, Hazan, and Ma]{GarberHM15}
	D.~Garber, E.~Hazan, and T.~Ma, ``Online learning of eigenvectors,'' in
	\emph{Proceedings of the 32nd International Conference on Machine Learning,
		{ICML} 2015, Lille, France, 6-11 July 2015}, 2015, pp. 560--568.
	
\end{thebibliography}

\end{document}

prelims.tex
\subsection{Notation}\label{sec:notation}
Bold lowercase letters such as $\u,\v,\w$ are used to denote vectors
and bold uppercase letters such as $\A,\B,\C$ to denote matrices. For
symmetric matrices $\A$ and $\B$, $\A \preceq \B$  denotes the
condition that $\x^\top \A \x \leq \x^\top \B \x$ for all $\x$ and
define $\B \succeq \A$ analogously. A symmetric matrix $\A$
is positive semidefinite if $\A \succeq \mZero$. For symmetric matrices
$\A, \B$, define their inner product as
$\iprod{\A}{\B} \eqdef \trace{\A^\top \B}$.

\subsection{Paper Organization}
The rest of this paper is organized as
follows. Section~\ref{sec:preliminaries} introduces  basic
mathematical facts used throughout the paper and also provides a proof
of the error bound of the standard batch method
(Theorem~\ref{thm:bernstein_kahan}). Section~\ref{sec:approach}
provides an overview of our approach to analyzing Oja's algorithm and
provides the main technical result of the paper. This technical result
is used in Section~\ref{sec:results} to prove the running time for
Oja's algorithm and to justify the choice of step size.
%, and further elucidate the performance improvements our algorithm achieves in different parameter regimes. 
Section~\ref{sec:proofs} presents the proof of the main technical
result. Section~\ref{sec:conc} concludes  and mentions a few interesting future directions. 

\section{Preliminaries}
\label{sec:preliminaries}


%%% Local Variables: 
%%% mode: latex
%%% TeX-master: "Ojasalgo"
%%% End: 
%
%\section{Setting}\label{sec:prelim}
%We receive (not necessarily symmetric)  matrices $\At{1},\At{2},\cdots,\At{n} \in \R^{d\times d}$ sampled independently satisfying the following properties for all $t \in [n]$:
%\begin{enumerate}
%	\item	$\E{\At{t}}=\Cov $ where $\Cov$ is symmetric
%	\item	$\twonorm{\At{t}- \Cov} \leq M$ and $\twonorm{\At{t}^\top - \Cov} \leq M$ with probability $1$, and
%	\item	$\max\left\{ \twonorm{\E{\left(\At{t}-\Cov\right) \trans{\left(\At{t}-\Cov\right)}}}, \twonorm{\E{ \trans{\left(\At{t}-\Cov\right)}\left(\At{t}-\Cov\right)}} \right\} \leq V$.
%\end{enumerate}
%We let $\v_1, ..., \v_d$ denote the eigenvectors of $\Cov$ and we let $\lambda_1 > \lambda_2 \geq \lambda_3 \geq ... \geq \lambda_d$ denote the corresponding eigenvalues. 
%
%Our goal is to estimate $\v_1$ given access to $\At{1},\cdots,\At{n}$. We do so by Oja's algorithm given in Algorithm~\ref{algo:Oja}.
%\begin{algorithm}[t]
%	\caption{Oja's algorithm for computing top eigenvector}\label{algo:Oja}
%	\begin{algorithmic}
%		\renewcommand{\algorithmicrequire}{\textbf{Input: }}
%		\renewcommand{\algorithmicensure}{\textbf{Output: }}
%		\REQUIRE $\At{1},\cdots,\At{n}$, initial vector $\wt{0}$.
%		\ENSURE  $\wt{n}$
%		\FOR{$t = 1,\cdots,n$}
%		\STATE $\wt{t} \leftarrow \wt{t-1} + \eta_t \At{t}\wt{t-1}$
%		\STATE $\wt{t} \leftarrow \wt{t}/\twonorm{\wt{t}}$
%		\ENDFOR
%		\STATE \textbf{Return} $\wt{n}$.
%	\end{algorithmic}
%	\label{algo:oja}
%\end{algorithm}
%
%%\textbf{Potential function}: Given a vector $\w$, we measure its distance from $\ustar$ using a potential function which is the $\sin$ of the angle between $\w$ and $\ustar$.
%%\begin{align*}
%%	\phiustar{\w} \eqdef \frac{\twonorm{\w}^2 - \iprod{\ustar}{\w}^2}{\twonorm{\w}^2}.
%%\end{align*}

The following basic inequalities regarding power
series, the exponential, and PSD matrices are used throughout. The
facts are summarized here:

\begin{lemma}[Basic Inequalities]\label{lem:trace-cauchy-schwarz}
The following are true:
\begin{itemize}
	\item	$1+x \leq \exp(x)$ for all $x$
	\item	$1+x \geq \exp\left(x-x^2\right)$ for all $x \geq 0$
	\item	$\frac{1}{1 + x} \leq \sum_{i=1}^{\infty} \frac{1}{(x+i)^2} \leq \frac{1}{x}$
	\item   $\iprod{\A}{\B} \leq \iprod{\A}{\C}$ for PSD matrices $\A, \B, \C$ with $\B \preceq \C$ 
	\item	$\trace{\A^\top \B} \leq \frac{1}{2} \trace{\A^\top \A + \B^\top \B}$ for all matrices $\A,\B \in \R^{m \times n}$.
\end{itemize}

\begin{proof}
The first inequality follows from the Taylor expansion of $\exp(x)$. The second comes from $1 + 0 = \exp(0 - 0^2)$ and $\frac{d}{dx} (1 + x) \leq \frac{d}{dx} \exp(x - x^2)$ for $x \geq 0$. The third follows by considering upper and lower Riemann sums of $\int_{y=1}^{\infty} 1/(x+y)$.
%simplifying: $\sum_{i=1}^{\infty} \int_{y=i}^{y=i+1}\frac{1}{(x+y)^2}dy\leq \sum_{i=1}^{\infty} \frac{1}{(x+i)^2}\leq \sum_{i=1}^{\infty} \int_{y=i-1}^{y=i}\frac{1}{(x+y)^2}dy$. 
The fourth from the fact that since $\A$ is PSD there is a matrix $\D$ with $\D^\top \D = \A$ and therefore 
\[
\iprod{\A}{\B} = \trace{\A^\top \B} = \trace{\D \B \D^\top}
\leq \trace{\D \C \D^\top} = \iprod{\A}{\C} ~.
\]
The final follows from Cauchy Schwarz and Young's inequality, i.e. $x\cdot y \leq \frac{1}{2}(x^2 + y^2)$ as
\[
\trace{\B^\top \A}
= \sum_{i \in [n]} 1_i \B^\top \A 1_i
\leq \sum_{i \in [n]} \twonorm{\A 1_i} \cdot \twonorm{\B 1_i} 
\leq \frac{1}{2} \sum_{i \in [n]} \left(\twonorm{\A 1_i}^2 + \twonorm{\B 1_i}^2\right)
\]
\end{proof}

\end{lemma}

%\subsection{Sample Complexity By Bernstein and Davis-Kahan}\label{sec:mb}
The following is a matrix Bernstein based proof of the error bound of the batch method. 
\begin{proof}[Proof of Theorem~\ref{thm:bernstein_kahan}]
Using Theorem 1.4 of \cite{Tropp12}, we have (w.p. $\geq 1-\delta$): 
\begin{equation}
\twonorm{\frac{1}{n}\sum_{i=1}^n \At{i} - \Cov } \le 2 \cdot \max\left\{\sqrt{\frac{\variance}{n}\log \frac{d}{\delta}}, \frac{\upperbound}{n}\log \frac{d}{\delta}\right\}.\label{eq:mb_1}
\end{equation}
Let $\widehat{\v}$ be the top eigenvector of $\widehat{\Cov}=\frac{1}{n}\sum_{i=1}^n \At{i}$. Using Wedin's theorem \cite{wedin1972perturbation}, implies:
\begin{equation}
\sin^2\langle\v_1, \widehat{\v}\rangle
\le \frac{\twonorm{\frac{1}{n}\sum_{i=1}^n \At{i} - \Cov}^2}{|\lambda_1-\lambda_2|^2 }. \label{eq:mb_2}
\end{equation}
Theorem now follows by combining \eqref{eq:mb_1} and \eqref{eq:mb_2}. \end{proof}


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

proof_thm.tex
%\begin{proof}[Proof of Theorem~\ref{thm:main}]
%The proof is essentially to substitute our choice of parameters in the result of Theorem~\ref{thm:main1}. Recall that Theorem~\ref{thm:main1} gives a bound of
%\begin{align}
%	\frac{1}{\delta\cdot Q} \cdot \left(d \cdot \exp\left(-2(\lambda_1-\lambda_2)\sum_{t \in [n]} \eta_t  + 5 \overline{\variance}\sum_{t\in[n]}\eta_t^2\right) + \variance \sum_{t=1}^{n} \eta_t^2 \exp
%	(-\sum_{j =t+1 }^n 2\eta_j(\lambda_1-\lambda_2) )\right), \label{eqn:master}
%\end{align}
%where $Q= \left( 1 - \frac{1}{\sqrt{\delta}}\sqrt{ \exp\left(18\overline{\variance}\sum_{t=1}^{n} \eta_i^2 \right) -  1}\right)$.
%%Recall our choice of stepsize
%%	\begin{align*}
%%	\eta_t = \left\{
%%	\begin{array}{l}
%%	\frac{\alpha}{(\lambda_1-\lambda_2)\beta \log d} \mbox{ for } t \leq \beta \log^2 d	 \\
%%	\frac{\alpha}{(\lambda_1-\lambda_2)(\beta+t)} \mbox{ for } t > \beta \log^2 d
%%	\end{array}
%%	\right..
%%	\end{align*}
%First, we bound the sum of $\eta_i^2$ as follows,
%\begin{align*}
%	\variancetotal \sum_{t\in [n]} \eta_t^2 = \variancetotal \left(\sum_{t\in[\beta \log^2 d
%		]} \eta_t^2 + \sum_{t > \beta \log^2 d}^{n} \eta_t^2\right)
%	\leq \frac{\alpha^2 \variancetotal}{(\lambda_1-\lambda_2)^2\beta} + \frac{\alpha^2\variancetotal}{(\lambda_1-\lambda_2)^2 \beta \log^2 d} \leq \frac{1}{20} \log(1+\delta/4).
%\end{align*}
%This means that $\exp\left(18 \variancetotal \sum_{t\in[n]}\eta_t^2\right) \leq 1+\frac{\delta}{4}$. So we have $Q \geq \frac{1}{2}$. On the other hand, we also have the following bound on the sum of $\eta_t$
%\begin{align}
%	2(\lambda_1-\lambda_2) \sum_{t\in [\beta \log^2 d]} \eta_t
%	\geq 2\alpha \log d \quad 
%	\mbox{ and } \quad 
%	2(\lambda_1-\lambda_2) \sum_{t > \beta \log^2 d}^n \eta_t
%	\geq 2\alpha \log \left(\frac{n}{\beta\log^2 d}\right).
%	\label{eqn:bound_sum_eta}
%\end{align}
%This gives us the following bound on the first term on the right hand side of~\eqref{eqn:master}.
%\begin{align}
%	d \cdot \exp\left(-2(\lambda_1-\lambda_2)\sum_{t \in [n]} \eta_t  + 5 \overline{\variance}\sum_{t\in[n]}\eta_t^2\right) &\leq 2 d \exp\left(-2\alpha \left(\log d + \log \left(\frac{n}{\beta \log^2 d}\right)\right)\right) \leq \left(\frac{\beta}{n}\right)^{2\alpha}. \label{eqn:plug1}
%\end{align}
%For the second term on the right hand side of~\eqref{eqn:master}, we break the summation into two parts
%\begin{align*}
%	\variance \sum_{t=1}^{n} \eta_t^2 \exp
%	(-\sum_{j =t+1 }^n 2\eta_j(\lambda_1-\lambda_2) )
%	&= \variance \sum_{t=1}^{\beta \log^2 d} \eta_t^2 \exp
%	(-\sum_{j =t+1 }^n 2\eta_j(\lambda_1-\lambda_2) ) \\
%	&\quad + \variance \sum_{t > \beta \log^2 d}^{n} \eta_t^2 \exp
%	(-\sum_{j =t+1 }^n 2\eta_j(\lambda_1-\lambda_2) ).
%\end{align*}
%We now bound each of the above terms.
%\begin{align}
%	&\variance \sum_{t=1}^{\beta \log^2 d} \eta_i^2 \exp
%	(-\sum_{j =t+1 }^n 2\eta_j(\lambda_1-\lambda_2) ) \nonumber\\
%	&= \variance \cdot \left(\frac{\alpha}{(\lambda_1-\lambda_2)\beta \log d}\right)^2 \exp
%	\left(-\sum_{j > \beta \log^2 d }^n 2\eta_j(\lambda_1-\lambda_2) \right) \sum_{t\in[\beta \log^2 d]} \exp\left(\frac{-2\alpha t}{\beta \log d}\right) \nonumber\\
%	&\leq \variance \left(\frac{\alpha}{(\lambda_1-\lambda_2)\beta \log d}\right)^2 \left(\frac{\beta\log^2 d}{n}\right)^{2\alpha} \cdot \frac{1}{1-\exp\left(-2\alpha/\beta \log d\right)} \quad (\mbox{ by~\eqref{eqn:bound_sum_eta}}) \nonumber\\
%	&\leq \variance \left(\frac{\alpha}{(\lambda_1-\lambda_2)\beta \log d}\right)^2 \left(\frac{\beta\log^2 d}{n}\right)^{2\alpha} \cdot \frac{\beta \log d}{\alpha} \quad (\mbox{ since } \frac{\alpha}{\beta \log d} < \frac{1}{4}) \nonumber\\
%	&\leq \left(\frac{\beta}{n}\right)^{2\alpha} \left(\log d\right)^{4\alpha -1}.\label{eqn:plug2}
%\end{align}
%On the other hand, we have
%\begin{align}
%	&\variance \sum_{t > \beta \log^2 d}^{n} \eta_t^2 \exp
%	(-\sum_{j =t+1 }^n 2\eta_j(\lambda_1-\lambda_2) ) \nonumber \\
%	&= \frac{\alpha^2 \variance}{(\lambda_1-\lambda_2)^2} \sum_{t > \beta \log^2 d}^{n} \frac{1}{t^2} \exp
%	(-2 \alpha \sum_{j =t+1 }^n \frac{1}{j} ) \nonumber \\
%	&\leq \frac{\alpha^2 \variance}{(\lambda_1-\lambda_2)^2} \sum_{t > \beta \log^2 d}^{n} \frac{1}{t^2} \exp
%	(-2 \alpha \log \frac{n}{t}) \nonumber \\
%	&= \frac{\alpha^2 \variance}{(\lambda_1-\lambda_2)^2} \sum_{t > \beta \log^2 d}^{n} \frac{1}{t^2} \left(\frac{t}{n}\right)^{2 \alpha}
%	\leq \frac{\alpha^2 \variance}{(\lambda_1-\lambda_2)^2} \cdot \frac{1}{n^{2\alpha}}\sum_{t > \beta \log^2 d}^{n} 
%	t^{2\alpha-2} \nonumber \\
%	&\leq \frac{2 \alpha \variance}{(\lambda_1-\lambda_2)^2 n} \quad (\mbox{ since } \alpha > 1/2 \mbox{ and } \sum_{x=1}^n x^\gamma \leq x^{\gamma+1}/(\gamma+1) \mbox{ for } \gamma > -1).
%	\label{eqn:plug3}
%\end{align}
%Plugging in~\eqref{eqn:plug1},~\eqref{eqn:plug2} and~\eqref{eqn:plug3} into~\eqref{eqn:master} proves the theorem.
%\end{proof}

\begin{proof}
	Recall that Theorem~\ref{thm:main1} gives a bound of
	\begin{align}\label{eqn:master}
		\frac{1}{Q} \exp\left(5 \overline{\variance}\sum_{i\in[n]}\eta_i^2\right)
		\left(d \cdot \exp\left(-2(\lambda_1-\lambda_2)\sum_{i \in [n]} \eta_i  \right) + \variance \sum_{i=1}^{n} \eta_i^2 \exp
		\left(-\sum_{j =i+1 }^n 2\eta_j(\lambda_1-\lambda_2) \right)\right)
%		\frac{d \cdot \gamma}{Q} \cdot \exp\left(\sum_{t \in [n]} (-2(\lambda_1-\lambda_2) \eta_t  + 5 
%		\variancetotal \eta_t^2)\right) + 
%		\frac{\variance \cdot \gamma}{Q} \sum_{t=1}^{n} \eta_t^2 \exp
%		\left(-\sum_{j =t+1 }^n 2\eta_j(\lambda_1-\lambda_2) \right), 
	\end{align}
	where $Q \eqdef \frac{\delta^2}{C\log(1/\delta)}\left( 1 - \frac{1}{\sqrt{\delta}}\sqrt{ \exp\left(18\overline{\variance}\sum_{i=1}^{n} \eta_i^2 \right) -  1}\right)$. % and $\gamma = C/(\delta^2 \log(1/\delta))$ for some constant $C$.	
	Since $\eta_i=\frac{\alpha}{(\lambda_1-\lambda_2)(\beta + i)}$, we have $\sum_{i\in [n]} \eta_i^2 \leq \frac{\alpha^2}{(\lambda_1-\lambda_2)^2\beta}$ and by our assumption that $\frac{\overline{\variance} \alpha^2}{(\lambda_1-\lambda_2)^2\beta} \leq \frac{1}{18}\log \left(1+\frac{\delta}{100}\right)$, we have: 
	\begin{equation}
	\label{eqn:main4}
	\exp\left(18 \overline{\variance}\sum_{i\in[n]}\eta_i^2\right) \leq \sqrt{2} \quad \Rightarrow \quad 
	Q\geq \frac{\delta^2}{C \log(1/\delta)}.
	\end{equation}
	Moreover, since $\sum_{i \in [n]} \eta_i \geq \frac{\alpha}{\lambda_1-\lambda_2} \log\left(1+n/\beta\right)$, we have
	\begin{align}
		\exp\left(-2(\lambda_1-\lambda_2)\sum_{i \in [n]} \eta_i\right) \leq \left(\frac{\beta}{\beta+n}\right)^{2\alpha}. \label{eqn:main4-5}
	\end{align}
	Note that $\sum_{j=i+1}^n \eta_j\leq \frac{\alpha}{\lambda_1-\lambda_2} \log\frac{n+\beta +1}{i+\beta+1}$. Moreover, as $\alpha > 1/2$, we have: 
	\begin{align}
	  &\sum_{i=1}^{n}\eta_i^2 \exp\left(-2(\lambda_1-\lambda_2)\sum_{j=i+1}^n \eta_j\right)\nonumber\\
	  &\leq \frac{\alpha^2}{(\lambda_1-\lambda_2)^2}\sum_{i=1}^{n} \frac{1}{(\beta+i)^2}\exp\left(2\alpha\log\frac{i+\beta+1}{n+\beta+1}\right),\notag\\&\leq \frac{(\beta+1)^2}{\beta^2}\cdot \frac{\alpha^2}{(\lambda_1-\lambda_2)^2(n+\beta+1)^{2\alpha}}\cdot \sum_{i=1}^n (i+\beta+1)^{2\alpha-2},\notag\\
	&\leq \frac{2\alpha^2}{(2\alpha-1)(\lambda_1-\lambda_2)^2(n+\beta+1)} \quad (\mbox{since } \alpha > 1/2 \mbox{ and } \sum_{i=1}^{n} i^{\gamma} \leq n^{\gamma+1}/(\gamma+1) \; \forall \; \gamma > -1). \label{eqn:main5}
	\end{align}
	Substituting~\eqref{eqn:main4},~\eqref{eqn:main4-5} and \eqref{eqn:main5} into~\eqref{eqn:master} proves the theorem.
\end{proof}
related.tex
\subsection{Comparison with Existing Results}\label{sec:comp}
Here we compare our sample complexity bounds with existing analyses of various methods. Recall that the error of the estimate $\w$ is $\sin^2(\w, \v_1)=1-(\w^\top \v_1)^2$.

We consider three popular methods used for computing $\v_1$. The first one is the batch method which computes largest eigenvector of empirical covariance and uses Wedin's theorem with matrix Bernstein inequality (cf. Theorem~\ref{thm:bernstein_kahan}). The second method is Alecton, which is very similar to Oja's algorithm \cite{SaOR15}. Finally, consider a block-power method (BPM) \cite{HardtP14, MitgliakisCJ13} which divides samples into different blocks and applies power iteration to the empirical estimate from each block. See Table~\ref{tab:comp} for the comparison.

We stress that some of the results we compare to make different assumptions than Definition~\ref{def:streaming-PCA}. The bounds stated for them are our best attempt to adapt their bounds in the setting of Definition~\ref{def:streaming-PCA} (which is quite standard). The next paragraph provides a simple example, which demonstrates the improvement in our result as compared to existing work.
%Table~\ref{tab:comp} compares sample complexity bounds for various methods. All the bounds hold with at least a constant probability; we ignore all the constants for brevity. The bounds are obtained using corresponding results applied to the standard Bernstein setting (see Section~\ref{sec:}). 
%
%Note that the second term in our sample complexity bound can be made smaller by selecting a larger $p$. Moreover, our first term is $O(\log d)$ better than the standard Matrix Bernstein bound \cite{Tropp12}. Hence, for small $\epsilon$, we can obtain a slightly better bound than even the best known Matrix Bernstein result. Furthermore, our bounds are $O(d\log(\frac{d}{\epsilon}))$ better than that of Alecton method \cite{SaOR15} and $O(\frac{\lambda_1}{\Delta})\log(\frac{d}{\epsilon})$ better than the block power method analysis of \cite{HardtP14}. 

Let $\At{i}=\x_i \x_i^{\top}$, where $\x_i\in \R^d$ and $\x_i=\e_1$ with probability $1/d$ and $\x_i=\sigma \e_j, 1\leq j\leq d$ with probability $1/d$ where $\e_j$ denotes the $j^{\textrm{th}}$ standard basis vector and $\sigma<1$. Note that $\Cov=\E{\At{i}}=\frac{(1-\sigma^2)}{d}\e_1\e_1^T+\frac{1}{d}\sigma^2 \eye$, $\|\At{i}\|_2\leq 1$ for all $i$, and $\|\E{\At{i}\At{i}^{\top}}\|_2\leq \frac{1}{d}$. Even for constant accuracy $\epsilon = \Omega(1)$, Theorem~\ref{thm:intro1} tells us that $n = \order{\frac{d \log^2 d}{(1-\sigma^2)^2}}$ is sufficient. On the other hand, Theorem $1$ of \cite{SaOR15} requires $n = \order{\frac{d^2 \log^2 d}{(1-\sigma^2)^2}}$, while Theorem $2.4$ of \cite{HardtP14} requires $n = \order{\frac{d \log^2 d}{(1-\sigma^2)^3}}$. Asymptotically, as $n$ becomes larger, our error scales as $\order{\frac{d}{(1-\sigma^2)^2 }\cdot \frac{1}{n}}$ while that of~\cite{SaOR15} scales as $\order{\frac{d^2}{(1-\sigma^2)^2 }\cdot \frac{\log n}{n}}$ and that of~\cite{HardtP14} scales as $\order{\frac{d}{(1-\sigma^2)^3 }\cdot \frac{\log n}{n}}$. Combining matrix Bernstein and Wedin's theorems gives an asymptotic error of $\order{\frac{d \log d}{(1-\sigma^2)^2 }\cdot \frac{1}{n}}$.
%
%bounds become strictly better than that of both \cite{HardtP14} as well as \cite{SaOR15}; in fact we obtain even tighter bound than the standard matrix Bernstein based bound. 
%Now, using Theorem~\ref{thm:main} with $p=8$, we require $n\geq  \frac{d}{(1-\sigma^2)^2\cdot \epsilon}+\frac{d\log^2 d}{(1-\sigma^2)\cdot \epsilon^{1/4}}$. In contrast, Theorem 1 of \cite{SaOR15} requires $n\geq  \frac{d^2\log(d/\epsilon)}{(1-\sigma^2)^2\cdot \epsilon}$, and Theorem 2.4 of \cite{HardtP14} requires $n\geq \frac{d\log(d/\epsilon)}{(1-\sigma^2)^3\cdot \epsilon}$.% Hence, approximately a larger multiplicative factor of $d$ and the gap $\lambda_1/\Delta$, respectively from our result. 
\begin{table}
\centering
\begin{tabular}[h]{|C{5cm}|c|c|}
\hline
Algorithm & Error & $\order{d}$ space? \\
\hline
Oja's \textbf{(this work, Theorem~\ref{thm:main})} & $\order{\frac{\variance }{(\lambda_1-\lambda_2)^2} \cdot \frac{1}{n}}$ & Yes \\
\hline
Matrix Bernstein + Wedin's theorem (Theorem~\ref{thm:bernstein_kahan}) & $\order{\frac{\variance\log d}{(\lambda_1-\lambda_2)^2}\cdot \frac{1}{n}}$ %+\frac{M^2\log^2 d}{\Delta^2n^2}$
& No \\
\hline
Alecton \cite{SaOR15} & $\order{\frac{\variance d}{(\lambda_1-\lambda_2)^2}\cdot \frac{\log n}{n}}$
%$\frac{d\variance}{\Delta^2\cdot n}\cdot \log\left(\frac{d\cdot \Delta\cdot n}{\variance}\right)$
& Yes \\
\hline
Block Power Method \cite{HardtP14} & $\order{\frac{\variance \lambda_1 \log d}{(\lambda_1-\lambda_2)^3} \cdot \frac{\log n}{n}}$
%$\left(\frac{\variance\log d}{\Delta^2\cdot n}+\frac{M^2\log^2 d}{\Delta^2n^2}\right)\cdot \frac{\lambda_1}{\Delta}\log(\frac{d\cdot \Delta\cdot n}{\variance})$
& Yes \\
\hline
%  &Bernstein&Alecton \cite{SaOR15}&Block PM \cite{HardtP14}\\[2pt]
%\hline Samples& $\frac{p^2\variance}{\Delta^2\epsilon}+\frac{M\lambda_1\log^2d}{\Delta^2\epsilon^{2/p}}$&$\frac{\variance\log d}{\Delta^2\epsilon}+\frac{M\log d}{\Delta\sqrt{\epsilon}}$ & $\frac{d\variance}{\Delta^2}\cdot \log\left(\frac{d}{\epsilon}\right)$&$\left(\frac{\variance\log d}{\Delta^2\epsilon}+\frac{M\log d}{\Delta\sqrt{\epsilon}}\right)\cdot \frac{\lambda_1}{\Delta}\log(\frac{d}{\epsilon})$\\[2pt]\hline
\end{tabular}
\caption{Asymptotic error guaranteed by various methods under assumptions of Definition~\ref{def:streaming-PCA} with at least constant probability, and ignoring constant factors. Recall that the error is defined as $\sin^2(\w, \v_1) = 1-(\w^\top \v_1)^2$. Our analysis provides the optimal $1/n$ error decay rate as compared to Alecton and Block power method which obtain $\frac{\log n}{n}$. Moreover, our bound is $O(d)$ tighter than that of Alecton \cite{SaOR15} and $O(\frac{\lambda_1}{\lambda_1-\lambda_2})$ tighter bound than that of Block Power Method \cite{HardtP14}. The assumptions made in \cite{SaOR15} for Alecton are different from our (more standard) assumption; we have optimized their bounds are optimized in our setting. See Section~\ref{sec:comp} for a concrete example where our analysis provides these improvements over \cite{SaOR15,HardtP14}.}\label{tab:comp}
\end{table}

\subsection{Additional Related Work}
Existing results for computing largest eigenvector of a data covariance matrix using streaming samples can be divided into three broad settings: a) stochastic data, b) arbitrary sequence of data, c) regret bounds for arbitrary sequence of data.

{\bf Stochastic data}: Here, the data is assumed to be sampled i.i.d. from a fixed distribution. The analysis of Oja's algorithm as well as those of block power method and Alecton mentioned earlier are in this setting. \cite{MitgliakisCJ13} also obtained a result in the restricted spiked covariance model.
%Moreover, their method is block power method and requires a multiplicative factor of $O(\lambda_1/\Delta)$ more samples than our result.
\cite{BalsubramaniDF13}  provides an analysis of a modification of Oja's algorithm but with an extra $O(d^5)$ multiplicative factor compared to ours.~\cite{jin2015robust} provides an algorithm based on shift and invert framework that obtains the same asymptotic error as ours. However, their algorithm requires warm start with a vector that is already constant close to the top eigenvector, which itself is a hard problem.

{\bf Arbitrary data}: In this setting, each row of the data matrix is provided in an arbitrary order. Most of the existing methods here first compute a sketch of the matrix and use that to compute an estimate of the top eigenvector \cite{ClarksonW09, Liberty13, nelson2013osnap, cohen2015optimal, GhashamiLPW15, BoutsidisGKL15}. However, a direct application of such techniques to the stochastic setting leads to sample complexity bounds which are larger by a multiplicative factor of $O(d)$ (ignoring other factors like variance etc). Finally, \cite{Shamir15, garber2015fast,jin2015robust} also provide methods for eigenvector computation, but they require multiple passes over the data and hence do not apply to the streaming setting. 


{\bf Regret bounds}: Here, at each step the algorithm has to output an estimate $\w$ of $\v_1$ for which we get reward of $\w^T \A_i \w$ and the goal is to minimize the regret w.r.t. $\v_1$. The algorithms in this regime are mostly based on online convex optimization and applying them in our setting would again result in a loss of multiplicative $O(d)$. Moreover, typical algorithms in this setting are not memory efficient \cite{WarmuthK06, GarberHM15}. 


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

results.tex
\section{Main Results}\label{sec:results}

Theorem~\ref{thm:main1}, from the previous section, leads to our main results, provided here. The theorem and proof are below and essentially consist of choosing appropriate parameters to efficiently apply Theorem~\ref{thm:main1}. Once we have this theorem, Theorems~\ref{thm:intro1} and~\ref{thm:intro2} follow by choosing $\alpha = \log d$ and $\alpha = 6$ respectively.


%\begin{theorem}\label{thm:main}
%	Fix any $\delta > 0$. Suppose the stepsize is chosen to be
%	\begin{align*}
%	\eta_t = \left\{
%		\begin{array}{l}
%			\frac{\alpha}{(\lambda_1-\lambda_2)\beta \log d} \mbox{ for } t \leq \beta \log^2 d	 \\
%			\frac{\alpha}{(\lambda_1-\lambda_2)t} \mbox{ for } t > \beta \log^2 d,
%		\end{array}
%		\right.
%	\end{align*}
%	where $\alpha > \frac{1}{2}$ and
%	\begin{align*}
%		\beta \eqdef 20\max\left(\frac{ \upperbound \alpha}{(\lambda_1-\lambda_2)},\frac{\left(\variance+\left(\lambda_{1}\right)^2\right) \alpha^2}{(\lambda_1-\lambda_2)^2\log \left(1+\frac{\delta}{100}\right)}\right).
%	\end{align*}
%	Suppose the number of samples $n > \beta \log^2 d$.
%	Then the output $\wt{n}$ of Algorithm~\ref{algo:Oja} satisfies:
%	\begin{align*}
%		1-(\wt{n}^\top \v_1)^2 \leq \frac{2d}{\delta} \left(\frac{\beta}{n}\right)^{2 \alpha \left(\lambda_{1}-\lambda_{2}\right)} + \frac{2 \alpha \variance}{(\lambda_1-\lambda_2)\delta}\cdot \frac{1}{n},
%	\end{align*}
%	\begin{align*}
%		1-(\wt{n}^\top \v_1)^2 \leq \frac{2 \alpha \variance}{(\lambda_1-\lambda_2)^2 n} + 2 \left(\frac{\beta}{n}\right)^{2\alpha} \left(\log d\right)^{4\alpha -1},
%	\end{align*}
%	with probability greater than $1-\delta$.
%\end{theorem}
\begin{theorem}\label{thm:main}
	Fix any $\delta > 0$ and suppose the step sizes are set to 
	$\eta_t = \frac{\alpha}{(\lambda_1-\lambda_2)(\beta+t)}$ for $\alpha > \frac{1}{2}$ and
	\begin{align*}
	\beta \eqdef 20\max\left(\frac{ \upperbound \alpha}{(\lambda_1-\lambda_2)},\frac{\left(\variance+\left(\lambda_{1}\right)^2\right) \alpha^2}{(\lambda_1-\lambda_2)^2\log \left(1+\frac{\delta}{100}\right)}\right).
	\end{align*}
	Suppose the number of samples $n > \beta$.
	Then the output $\wt{n}$ of Algorithm~\ref{algo:oja} satisfies:
	\begin{align*}
	1-(\wt{n}^\top \v_1)^2 \leq \frac{C\log(1/\delta)}{\delta^2}\left({d} \left(\frac{\beta}{n}\right)^{2 \alpha} + \frac{\alpha^2 \variance}{(2\alpha-1)(\lambda_1-\lambda_2)^2}\cdot \frac{1}{n}\right),
	\end{align*}
	with probability at least $1-\delta$. Here $C$ is an absolute numerical constant.
\end{theorem}
\input{proof_thm}
%The error in above expression comprises of two parts: the first part can be thought of as error due to initial point and the second part can be thought of as asymptotic statistical error. Given a target accuracy $\epsilon$, we can use Theorem~\ref{thm:main} to choose appropriate parameter values and obtain the optimal number of samples. We will now provide two corollaries instantiating Theorem~\ref{thm:main} under two different parameter settings. Before we give those corollaries, we make the following notation indicating what we think of as samples required to overcome initial error and those required to achieve asymptotic statistical accuracy. Define
%%\begin{align*}
%%	\ninit &\eqdef \max\left(\frac{M}{\lambda_1 - \lambda_2}, \frac{\variance }{(\lambda_1-\lambda_2)^2}, \frac{\lambda_1^2 }{(\lambda_1-\lambda_2)^2}\right), \qquad
%%	\nfinal \eqdef \frac{\variance}{(\lambda_1-\lambda_2)^2}.
%%\end{align*}
%%give the following lemma which tells us the guarantee given by matrix Bernstein applied along with Davis-Kahan theorem.
%\praneeth{TODO: Put a version of matrix Bernstein here.}
%\begin{lemma}[Matrix Bernstein+Davis-Kahan]\label{lem:matrixBern}
%	\begin{align*}
%	n > c \left(\frac{\variance \log {d}}{(\lambda_1 -\lambda_2)^2} \cdot \frac{1}{\epsilon} + \frac{M\log {d}}{\lambda_1 -\lambda_2} \cdot \frac{1}{\sqrt{\epsilon}}\right).
%	\end{align*}
%\end{lemma}
%We now give our first result, which matches Lemma~\ref{lem:matrixBern} up to a multiplicative factor of $\log d$ and an extra additive term.
%\begin{corollary}\label{cor:1}
%	There exist absolute numerical constants $c_1$, $c_2$ and $c_3$ such that the following holds.
%	Fix any $\epsilon > 0$ and $p > 1$. Suppose the stepsize is chosen to be $\eta_t = \frac{A}{N+t}$, where $A = c_1 \frac{\log {d}}{\lambda_1 - \lambda_2}$, $N = c_2 \ninit \log^2 \frac{d}{\delta}$. If
%	\begin{align*}
%		n > c_3 \left(\frac{\variance \log^2 {d}}{ (\lambda_1-\lambda_2)^2} \cdot \frac{1}{\epsilon} + \max\left(\frac{M \log {d}}{\lambda_1 - \lambda_2}, \frac{\lambda_1^2 \log^2 {d}}{(\lambda_1-\lambda_2)^2}\right) \cdot \frac{1}{\sqrt{\epsilon}}\right),
%	\end{align*}
%	then after $n$ steps, with probability at least $\frac{9}{10}$, the error in Algorithm~\ref{algo:oja} is at most $1-(\wt{n}^\top \v_1)^2 \leq \epsilon$.
%\end{corollary}
%\praneeth{Point out somewhere that $\order{1/\epsilon}$ is higher order, rest lower order.}
%The additive lower order term of $\frac{\lambda_1^2 \log^2 {d}}{(\lambda_1-\lambda_2)^2}\cdot \frac{1}{\sqrt{\epsilon}}$ is due to computational reasons: even if the variance $\variance$ in our distribution were equal to zero i.e., we obtain the matrix deterministically in one sample, our streaming setup does not let us estimate the top eigenvector of the matrix in one step. Matrix Bernstein, on the other hand, being a purely statistical result does not consider the computational aspect and declares that we can estimate the top eigenvector using a single sample. It might be possible to reduce our bound on the additive lower order term but getting rid of it completely is not possible. Coming to the leading order terms with $\order{1/\epsilon}$, we note that our bound is off by a factor of $\log d$ as compared to Lemma~\ref{lem:matrixBern}. However, for small enough error $\epsilon$, we can do significantly better. More concretely, we have the following corollary of our main result.
%\begin{corollary}\label{cor:2}
%	Let $\epsilon < d^{-\alpha}$ for some $\alpha > 0$. There exist absolute numerical constants $c_1(\alpha)$, $c_2(\alpha)$ and $c_3(\alpha)$ such that the following holds.
%	Suppose the stepsize is chosen to be $\eta_t = \frac{A}{N+t}$
%	where $A = \frac{c_1}{\lambda_1 - \lambda_2}$, $N = c_2 \max(\frac{M}{\lambda_1 - \lambda_2}, \frac{\variance }{(\lambda_1-\lambda_2)^2}, \frac{\lambda_1^2 }{(\lambda_1-\lambda_2)^2})$. If	
%	\begin{align*}
%	n > c_3 \left( \frac{\variance}{(\lambda_1-\lambda_2)^2} \cdot \frac{1}{\epsilon} +
%	\max\left(\frac{M}{\lambda_1 - \lambda_2}, \frac{\lambda_1^2 }{(\lambda_1-\lambda_2)^2}\right) \cdot \frac{1}{\sqrt[10]{\epsilon}} \right),
%	\end{align*}
%	then after $n$ steps, with probability at least $\frac{9}{10}$, the error in Algorithm~\ref{algo:oja} is at most $1-(\wt{n}^\top \v_1)^2 \leq \epsilon$.
%\end{corollary}
%Note that in this low error regime, our result is significantly better than Lemma~\ref{lem:matrixBern}. Most strikingly, we improve upon the leading order term with $\order{1/\epsilon}$ in our result by a factor of $\order{\log d}$ as compared to Lemma~\ref{lem:matrixBern}. Moreover, our dependence on the maximum norm bound $M$ is better by a factor of $\frac{\log d}{\poly{\epsilon}}$. As discussed earlier, we still have an additive lower order term which arises due to computational reasons. %Comparing Corollary~\ref{cor:1} with Lemma~\ref{lem:matrixBern}, we see that in the large error regime, i.e., when $\epsilon = \Omega(1)$, we require the number of samples $n$ to grow as
%%\begin{align*}
%%	n > c \max\left(\frac{pM \log \frac{d}{\delta}}{\lambda_1 - \lambda_2}, \frac{p^2\variance \log^2 \frac{d}{\delta}}{(\lambda_1-\lambda_2)^2}, \frac{p^2\lambda_1^2 \log^2 \frac{d}{\delta}}{(\lambda_1-\lambda_2)^2}\right),
%%\end{align*}
%%which is larger than that required by Lemma~\ref{lem:matrixBern}. The extra term of $\frac{p^2\lambda_1^2 \log^2 \frac{d}{\delta}}{(\lambda_1-\lambda_2)^2}$ is due to the fact that our algorithm requires a new sample whenever it performs a computation, while the argument for matrix Bernstein uses the same samples for multiple steps of computation. Asymptotically, as $\epsilon \rightarrow 0$, our result is off from Bernstein by a factor of $\order{\log \frac{d}{\delta}}$. The following corollary, with a different setting of parameters fixes this issue.
%
%\begin{corollary}\label{cor:3}
%	Fix any $\delta > 0$, $\epsilon > 0$ such that $\epsilon < d^{-\alpha}$ for some $\alpha > 0$, and $p > \frac{2}{\alpha}$.
%	Suppose the stepsize is chosen to be $\eta_t = \frac{A}{N+t}$
%	where $A = \frac{p}{\lambda_1 - \lambda_2}$, $N = 40\max(\frac{pM}{\lambda_1 - \lambda_2}, \frac{p^2\variance }{(\lambda_1-\lambda_2)^2}, \frac{p^2\lambda_1^2 }{(\lambda_1-\lambda_2)^2})$. If	
%	\begin{align*}
%	n > \frac{4p^2 \variance}{(\lambda_1-\lambda_2)^2 \delta} \cdot \frac{1}{\epsilon} +
%	80 \max\left(\frac{pM}{\lambda_1 - \lambda_2}, \frac{p^2\variance }{(\lambda_1-\lambda_2)^2}, \frac{p^2\lambda_1^2 }{(\lambda_1-\lambda_2)^2}\right) \frac{1}{\delta^{1/p}} \cdot \left(\frac{1}{\epsilon}\right)^{2/p},
%	\end{align*}
%	then after $n$ steps, with probability at least $1-\delta$, the error in Algorithm~\ref{algo:oja} is at most
%	\begin{align*}
%	1-(\wt{n}^\top \v_1)^2 \leq \epsilon.
%	\end{align*}
%\end{corollary}
%Corollary~\ref{cor:2} shows that asymptotically, as $\epsilon \rightarrow 0$, the number of samples required by Oja's algorithm is \emph{better} by a factor of $\order{\log\frac{d}{\delta}}$ than the guarantee given by an offline application of matrix Bernstein and Davis-Kahan. This matches the recent results of~\cite{garber2015fast,jin2015robust}. However, note that these papers require an initial point with error smaller than a small constant, as compared to our result which does not require a good initial point.
%
%We would also like to point out that in both Corollaries~\ref{cor:1} and~\ref{cor:2}, the dependence of $n$ on $M$ is very mild. In particular, we only need $n > \otilde{\frac{M}{\epsilon^{1/p}}}$ for any $p > 1$, where as Lemma~\ref{lem:matrixBern} has a dependence of $n > \otilde{\frac{M}{\sqrt{\epsilon}}}$. It is easy to formulate settings where the sample complexity guaranteed by our results is significantly better than that guaranteed by Lemma~\ref{lem:matrixBern}. We do not fully know if the reason for this is a weakness in deriving Bernstein inequality or because it demands high probability.
%

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

1602.06929.extracted_bib
[Jolliffe(2002)]{jolliffe2002principal I.~Jolliffe, Principal component analysis.\hskip 1em plus 0.5em minus 0.4em\relax Wiley Online Library, 2002.
1602.06929.extracted_bib
[Hall et~al.(1998)Hall, Marshall, and Martin]{hall1998incremental P.~M. Hall, A.~D. Marshall, and R.~R. Martin, ``Incremental eigenanalysis for classification.'' in BMVC, vol.~98.\hskip 1em plus 0.5em minus 0.4em\relax Citeseer, 1998, pp. 286--295.
1602.06929.extracted_bib
[Weng et~al.(2003)Weng, Zhang, and Hwang]{weng2003candid J.~Weng, Y.~Zhang, and W.-S. Hwang, ``Candid covariance-free incremental principal component analysis,'' Pattern Analysis and Machine Intelligence, IEEE Transactions on, vol.~25, no.~8, pp. 1034--1040, 2003.
1602.06929.extracted_bib
[Ross et~al.(2008)Ross, Lim, Lin, and Yang]{ross2008incremental D.~A. Ross, J.~Lim, R.-S. Lin, and M.-H. Yang, ``Incremental learning for robust visual tracking,'' International Journal of Computer Vision, vol.~77, no. 1-3, pp. 125--141, 2008.
1602.06929.extracted_bib
[Krasulina(1970)]{krasulina1970method T.~Krasulina, ``Method of stochastic approximation in the determination of the largest eigenvalue of the mathematical expectation of random matrices,'' Automatation and Remote Control, pp. 50--56, 1970.
1602.06929.extracted_bib
[Oja(1982)]{oja1982simplified E.~Oja, ``Simplified neuron model as a principal component analyzer,'' Journal of mathematical biology, vol.~15, no.~3, pp. 267--273, 1982.
1602.06929.extracted_bib
[Balsubramani et~al.(2013)Balsubramani, Dasgupta, and Freund]{BalsubramaniDF13 A.~Balsubramani, S.~Dasgupta, and Y.~Freund, ``The fast convergence of incremental pca,'' in Advances in Neural Information Processing Systems, 2013, pp. 3174--3182.
1602.06929.extracted_bib
[Mitliagkas et~al.(2013)Mitliagkas, Caramanis, and Jain]{MitgliakisCJ13 I.~Mitliagkas, C.~Caramanis, and P.~Jain, ``Memory limited, streaming pca,'' in Advances in Neural Information Processing Systems, 2013, pp. 2886--2894.
1602.06929.extracted_bib
[Sa et~al.(2015)Sa, Re, and Olukotun]{SaOR15 C.~D. Sa, C.~Re, and K.~Olukotun, ``Global convergence of stochastic gradient descent for some non-convex matrix problems,'' in Proceedings of the 32nd International Conference on Machine Learning, {ICML 2015, Lille, France, 6-11 July 2015, 2015, pp. 2332--2341.
1602.06929.extracted_bib
[Vershynin(2010)]{vershynin2010introduction R.~Vershynin, ``Introduction to the non-asymptotic analysis of random matrices,'' arXiv preprint arXiv:1011.3027, 2010.
1602.06929.extracted_bib
[Tropp(2012)]{Tropp12 J.~A. Tropp, ``User-friendly tail bounds for sums of random matrices,'' Foundations of Computational Mathematics, vol.~12, no.~4, pp. 389--434, 2012.
1602.06929.extracted_bib
[Johnstone(2001)]{johnstone2001distribution I.~M. Johnstone, ``On the distribution of the largest eigenvalue in principal components analysis,'' Annals of statistics, pp. 295--327, 2001.
1602.06929.extracted_bib
[Wedin(1972)]{wedin1972perturbation P.-{\AA. Wedin, ``Perturbation bounds in connection with singular value decomposition,'' BIT Numerical Mathematics, vol.~12, no.~1, pp. 99--111, 1972.
1602.06929.extracted_bib
[Jin et~al.(2015)Jin, Kakade, Musco, Netrapalli, and Sidford]{jin2015robust C.~Jin, S.~M. Kakade, C.~Musco, P.~Netrapalli, and A.~Sidford, ``Robust shift-and-invert preconditioning: Faster and more sample efficient algorithms for eigenvector computation,'' arXiv preprint arXiv:1510.08896, 2015.
1602.06929.extracted_bib
[Golub and Van~Loan(2012)]{golub2012matrix G.~H. Golub and C.~F. Van~Loan, Matrix computations.\hskip 1em plus 0.5em minus 0.4em\relax JHU Press, 2012, vol.~3.
1602.06929.extracted_bib
[Hardt and Price(2014)]{HardtP14 M.~Hardt and E.~Price, ``The noisy power method: {A meta algorithm with applications,'' in Advances in Neural Information Processing Systems 27: Annual Conference on Neural Information Processing Systems 2014, December 8-13 2014, Montreal, Quebec, Canada, 2014, pp. 2861--2869.
1602.06929.extracted_bib
[Cohen et~al.(2010)Cohen, Lee, Miller, Pachocki, and Sidford]{cohen2016median M.~Cohen, Y.~T. Lee, G.~Miller, J.~Pachocki, and A.~Sidford, ``Geometric median in nearly linear time,'' To Appear in 48th Annual Symposium on the Theory of Computing (STOC) 2016, 2010.
1602.06929.extracted_bib
[Clarkson and Woodruff(2009)]{ClarksonW09 K.~L. Clarkson and D.~P. Woodruff, ``Numerical linear algebra in the streaming model,'' in Proceedings of the forty-first annual ACM symposium on Theory of computing.\hskip 1em plus 0.5em minus 0.4em\relax ACM, 2009, pp. 205--214.
1602.06929.extracted_bib
[Liberty(2013)]{Liberty13 E.~Liberty, ``Simple and deterministic matrix sketching,'' in Proceedings of the 19th ACM SIGKDD international conference on Knowledge discovery and data mining.\hskip 1em plus 0.5em minus 0.4em\relax ACM, 2013, pp. 581--588.
1602.06929.extracted_bib
[Nelson and Nguy{\^en(2013)]{nelson2013osnap J.~Nelson and H.~L. Nguy{\^en, ``Osnap: Faster numerical linear algebra algorithms via sparser subspace embeddings,'' in Foundations of Computer Science (FOCS), 2013 IEEE 54th Annual Symposium on.\hskip 1em plus 0.5em minus 0.4em\relax IEEE, 2013, pp. 117--126.
1602.06929.extracted_bib
[Cohen et~al.(2015)Cohen, Nelson, and Woodruff]{cohen2015optimal M.~B. Cohen, J.~Nelson, and D.~P. Woodruff, ``Optimal approximate matrix product in terms of stable rank,'' arXiv preprint arXiv:1507.02268, 2015.
1602.06929.extracted_bib
[Ghashami et~al.(2015)Ghashami, Liberty, Phillips, and Woodruff]{GhashamiLPW15 M.~Ghashami, E.~Liberty, J.~M. Phillips, and D.~P. Woodruff, ``Frequent directions: Simple and deterministic matrix sketching,'' arXiv preprint arXiv:1501.01711, 2015.
1602.06929.extracted_bib
[Boutsidis et~al.(2015)Boutsidis, Garber, Karnin, and Liberty]{BoutsidisGKL15 C.~Boutsidis, D.~Garber, Z.~Karnin, and E.~Liberty, ``Online principal components analysis,'' in Proceedings of the Twenty-Sixth Annual ACM-SIAM Symposium on Discrete Algorithms.\hskip 1em plus 0.5em minus 0.4em\relax SIAM, 2015, pp. 887--901.
1602.06929.extracted_bib
[Shamir(2015)]{Shamir15 O.~Shamir, ``A stochastic {PCA and {SVD algorithm with an exponential convergence rate,'' in Proceedings of the 32nd International Conference on Machine Learning, {ICML 2015, Lille, France, 6-11 July 2015, 2015, pp. 144--152.
1602.06929.extracted_bib
[Garber and Hazan(2015)]{garber2015fast D.~Garber and E.~Hazan, ``Fast and simple pca via convex optimization,'' arXiv preprint arXiv:1509.05647, 2015.
1602.06929.extracted_bib
[Warmuth and Kuzmin(2006)]{WarmuthK06 M.~K. Warmuth and D.~Kuzmin, ``Randomized {PCA algorithms with regret bounds that are logarithmic in the dimension,'' in Advances in Neural Information Processing Systems 19, Proceedings of the Twentieth Annual Conference on Neural Information Processing Systems, Vancouver, British Columbia, Canada, December 4-7, 2006, 2006, pp. 1481--1488.
1602.06929.extracted_bib
[Garber et~al.(2015)Garber, Hazan, and Ma]{GarberHM15 D.~Garber, E.~Hazan, and T.~Ma, ``Online learning of eigenvectors,'' in Proceedings of the 32nd International Conference on Machine Learning, {ICML 2015, Lille, France, 6-11 July 2015, 2015, pp. 560--568.