SNIPE for Memory-Limited PCA From Incomplete Data
Armin Eftekhari
Laura Balzano
Dehui Yang
Michael B. Wakin\footnote{AE is with the Alan Turing Institute in London. LB is with the Department of Electrical Engineering and Computer Science at the University of Michigan
Ann Arbor. DY and MBW are with the Department of Electrical Engineering and Computer Science at Colorado School of Mines. (E-mails: aeftekhari@turing.ac.uk; girasole@umich.edu; dyang@mines.edu; mwakin@mines.edu)
{The linear subspace model is pervasive in science and engineering and particularly in large datasets which are often incomplete due to missing measurements and privacy issues. Therefore, a critical problem in modeling is to develop algorithms for estimating a low-dimensional subspace model from incomplete data efficiently in terms of both computational complexity and memory storage. In this paper we study an algorithm that processes blocks of incomplete data to estimate the underlying subspace model. Our algorithm has a simple interpretation as optimizing the subspace to fit the observed data block but remain close to the previous estimate. We prove a linear rate of convergence for the algorithm and our rate holds with high probability.
SNIPE_v10.tex
 %% LyX 2.1.3 created this file.  For more info, see http://www.lyx.org/.
%% Do not edit unless you really know what you are doing.
\documentclass[english]{article}
\usepackage[T1]{fontenc}
\usepackage[latin9]{inputenc}
\usepackage{verbatim}
\usepackage{amsthm}
\usepackage{amsmath}
\usepackage{amssymb}
\usepackage{esint}
\usepackage{varwidth}
\usepackage{url}
\usepackage{constants}
\usepackage{graphicx}
\usepackage{subfig}
\usepackage{float}


\makeatletter
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Textclass specific LaTeX commands.
\theoremstyle{plain}
\newtheorem{thm}{\protect\theoremname}
  \theoremstyle{plain}
  \newtheorem{lem}[thm]{\protect\lemmaname}
  \theoremstyle{plain}
  \newtheorem{cor}[thm]{\protect\corollaryname}
 \theoremstyle{plain}
  \newtheorem{rem}{Remark}

% latin
\def\ie{{\em i.e.,~}}
\def\eg{{\em e.g.,~}}
\def\cf{{\em cf.,~}}
\def\ea{{\em et al.~}}

\newcommand{\note}[1]{\textbf{\textcolor{blue}{#1}}}
\newcommand{\laura}[1]{\textbf{\textcolor{magenta}{#1}}}

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% User specified LaTeX commands.
\usepackage{hyperref}
\def \SU {\mathcal{S}}
%\def \SU {\mathcal{S}}
\def \A {\mathcal{A}}
\def \B {\mathcal{B}}
\def \alg {\operatorname{SNIPE}}
\def \Q {\mathcal{Q}}
\def \E {\mathbb{E}}
\def \h {\widehat}
\def \U {\mathcal{U}}
\def \grouse {\operatorname{GROUSE}}
\def \sS {\mathcal{S}}

\hypersetup{
   colorlinks=true,%
   citecolor=black,%
   filecolor=black,%
   linkcolor=red,%
   urlcolor=blue
}


\def \l{\left}
\def \r{\right}

% 1 inch margins when printed from ps
\hoffset=0in \voffset=0in \evensidemargin=0in \oddsidemargin=0in
\textwidth=6.5in \topmargin=0in \headheight=0.0in \headsep=0.0in
\textheight=9in

\makeatother

\usepackage{babel}
  \providecommand{\corollaryname}{Corollary}
  \providecommand{\lemmaname}{Lemma}
\providecommand{\theoremname}{Theorem}

\begin{document}

\title{SNIPE
for Memory-Limited PCA From Incomplete Data
}
\author{Armin Eftekhari, Laura Balzano, Dehui Yang, Michael B. Wakin\footnote{AE is with the Alan Turing Institute in London. LB is with the Department of Electrical Engineering and Computer Science at the University of Michigan, Ann Arbor. DY and MBW are with the Department of  Electrical Engineering and Computer Science at Colorado School of Mines. (E-mails: aeftekhari@turing.ac.uk; girasole@umich.edu; dyang@mines.edu; mwakin@mines.edu)}}

\maketitle
\begin{abstract}

{The linear subspace model is pervasive in science and engineering and particularly  in large datasets which are often incomplete due to missing measurements and privacy issues. Therefore, a critical problem in modeling is to develop algorithms for estimating a low-dimensional subspace model from incomplete data efficiently in terms of both computational complexity and memory storage. In this paper we study an algorithm that processes blocks of incomplete data to estimate the underlying subspace model. Our algorithm has a simple interpretation as optimizing the subspace to fit the observed data block but remain close to the previous estimate. We prove  a linear rate of convergence for the algorithm and our rate holds with high probability.}

\end{abstract}

\section{Introduction}
\label{sec:problem statement}

Linear models are the backbone of modern sciences and \emph{principal component analysis} (PCA) has traditionally been an indispensable tool in studying  collected data \cite{van2012subspace,ardekani1999activation,krim1996two,tong1998multichannel}. In dimensionality reduction, for example, PCA searches for the linear model that best describes the data \cite{cunningham2015linear}.
In this work we are particularly interested in applying PCA to data that suffers from erasures, while limited storage is available. In recommender systems, for instance, data is highly incomplete and yet so massive that can only be processed in small chunks \cite{recht2011hogwild}.

To focus our efforts and as our guiding example throughout, consider incomplete data from an unknown subspace, presented sequentially to the user who, due to hardware limitations, can only store  small amounts of data. We are here interested in developing a \emph{streaming} algorithm for PCA from incomplete measurements.





More concretely, consider an $r$-dimensional subspace $\SU$ with orthonormal basis $S\in\mathbb{R}^{n\times r}$.
For an  integer $T$, let the coefficient vectors $\{q_t\}_{t=1}^T\subset \mathbb{R}^r$ be independent copies of a random vector $q\in\mathbb{R}^r$.
 At time $t\in[1:T]:=\{1,2,\cdots,T\}$,  we observe each entry of  $s_t := S\cdot q_t \in \SU$ independently with a probability of $p$, and we collect the measurements in $y_t\in\mathbb{R}^n$, supported on a random index set $\omega_t\subseteq [1:n]$. Formally, we will  write this measurement process as $y_t=P_{\omega_t}(s_t)=P_{\omega_t}\cdot s_t$, where  $P_{\omega_t} \in\mathbb{R}^{n\times n}$ is the projection onto the coordinate set $\omega_t$, i.e.\ it equals one on its diagonal entries corresponding to the index set $\omega_t$, and is zero elsewhere.


Our objective is to identify the subspace $\SU$ from the measurements $\{y_t\}_{t=1}^T$ supported on the index sets $\{\omega_t\}_{t=1}^T$. Assuming that $r=\mbox{dim}(\SU)$ is known \emph{a priori} (or estimated from data by other means), we present the $\alg$ algorithm for this task in Section \ref{sec:proposed alg} and, in particular, provide helpful insight about the algorithm in Section \ref{sec:interp}. As summarized in Section \ref{sec:theory}, $\alg$ converges, globally and linearly, to the true subspace under  reasonable requirements. A detailed review of prior art is deferred to Section~\ref{sec:related work}. The performance of $\alg$  and rival algorithms are examined numerically in Section \ref{sec:simulations}, where we find that $\alg$  improves upon the state of the art.  Technical proofs appear in Section \ref{sec:analysis} and in the appendices, with Appendix \ref{sec:Toolbox} (Toolbox) collecting some of the frequently used mathematical tools.

Before delving into the details, let us conclude this section with an example. Suppose that $\SU\subset\mathbb{R}^{1000}$ is a generic subspace of dimension $r=3$ and take  $p=0.1$. Then, $\alg$ produces a sequence of estimates of $\SU$; the estimation error (with a metric to be specified shortly) versus $t$ is plotted in Figure \ref{fig:exampleIntro}.
\begin{figure}[H]
\begin{center}
\includegraphics[width=0.49\textwidth]{example_intro.pdf}
\caption{\label{fig:exampleIntro}This paper introduces $\alg$ for memory-limited PCA from incomplete data. This figure shows the estimation error of $\alg$ versus time in recovering a generic $3$-dimensional subspace from data subsampled by a factor of $10$, received sequentially. The algorithm and the error metric are described below.}
\end{center}
\end{figure}














\section{$\alg$}
\label{sec:proposed alg}

In this section, we propose Subspace Navigation  via Interpolation from Partial Entries (SNIPE), an algorithm for subspace identification from incomplete data, received sequentially.

 To formally present $\alg$, let us introduce some additional notation.
Recall the incoming sequence of measurement vectors $\{y_t\}_{t=1}^T\subset\mathbb{R}^n$, supported on the index sets $\{\omega_t\}_{t=1}^T\subseteq [1:n]$.
For an integer $K$, fix block sizes $\{b_k\}_{k=1}^K$ such that $\sum_{k=1}^K b_k = T$ and $b_k\ge r$ for all $k$.  We partition the data into $K$ (non-overlapping) blocks $\{Y_k\}_{k=1}^K$, where $Y_k\in\mathbb{R}^{n\times b_k}$ for each $k$.



 At a high level, $\alg$ processes the first measurement block  $Y_1$ to produce an estimate $\h{\SU}_1$ of the true subspace $\SU$. This estimate is then iteratively updated  after receiving each of the new blocks $\{Y_k\}_{k= 2}^K$,
 thereby producing a sequence of estimates $\{\h{\SU}_k\}_{k= 2}^K$. Every $\h{\SU}_k$ is an $r$-dimensional subspace of $\mathbb{R}^n$ with orthonormal basis $\h{S}_k\in\mathbb{R}^{n\times r}$; the particular choice of orthonormal basis is inconsequential throughout the paper.


More concretely, $\alg$ sets $\h{\SU}_1$ to be the span of the top $r$ left singular vectors of $Y_1$, namely the left singular vectors corresponding to the largest $r$ singular values of $Y_1$, with ties broken arbitrarily. Then, at iteration $k\in[2:K]$ and given the previous estimate $\h{\SU}_{k-1}=\mbox{span}( \h{S}_{k-1})$,  $\alg$ processes the columns of the $k$th measurement block $Y_k$ by forming the matrix
\begin{equation}
R_k =
\l[
\begin{array}{ccc}
\cdots
&
y_t + P_{\omega_t^C} \h{S}_{k-1}\l( P_{\omega_t} \h{S}_{k-1} \r)^{\dagger} y_t
&
\cdots
\end{array}
\r]
\in\mathbb{R}^{n\times b_k},
\qquad  \sum_{k'=1}^{k-1} b_{k'}+1 \le t \le \sum_{k'=1}^k b_{k'},
\label{eq:Rk}
\end{equation}
where $P_{\omega_t^C}=I_n-P_{\omega_t}\in\mathbb{R}^{n\times n}$ projects a vector onto the complement of the index set $\omega_t$. $\alg$ then updates its  estimate by setting $\h{\SU}_k$ to be the span of the top $r$ left singular vectors of $R_k$. Figure \ref{fig:Alg} summarizes these steps.

\begin{figure}
\begin{center}%
\fbox{\begin{minipage}[t]{1\columnwidth}%
\textbf{Input:}
\begin{itemize}
\item Dimension $r$,
\item Received data $\{y_t\}_{t=1}^T\subset \mathbb{R}^n$ supported on index sets $\{\omega_t\}_{t=1}^T\subseteq [1:n]$,
\item  Number of blocks $K$ and block sizes $\{b_k\}_{k=1}^K$, where $\sum_{k=1}^K b_k=T$ and $b_k\ge r$ for all $k$.
\end{itemize}

\textbf{Output: } $r$-dimensional subspace $\widehat{\SU}_K$.\\


\textbf{Body:}
\begin{itemize}


\item Form $Y_1\in\mathbb{R}^{n\times b_1}$ from the first $b_1$ measurement vectors $\{y_t\}_{t=1}^{b_1}$. Let $\h{\SU}_1$, with orthonormal basis $\h{S}_1\in\mathbb{R}^{n\times r} $,  be the span of the top $r$ left singular vectors of $Y_1$, namely those corresponding to the largest singular values. Ties are broken arbitrarily.
\item For $k \in [2:K]$, repeat:
\begin{itemize}
\item Set $R_{k}\leftarrow \{\}$.
\item For $t= \sum_{k'=1}^{k-1}b_{k'}+1:\sum_{k'=1}^{k} b_{k'}$, repeat
\begin{itemize}
\item  Set $$R_{k}\leftarrow \l[
\begin{array}{cc}
R_{k}
&
y_{t} +
P_{\omega_{t}^C} \h{S}_{k-1}  \l( P_{\omega_t} \h{S}_{k-1}\r)^{\dagger}  y_{t}
\end{array}
\r],
$$
where  $P_{\omega_t} \in\mathbb{R}^{n\times n}$ equals one on its  diagonal entries corresponding to the index set $\omega_t$, and is zero elsewhere. Likewise, $P_{\omega_t^C}$ projects a vector onto the complement of the index set $\omega_t$.
\end{itemize}
\item Let $\h{\SU}_{k}$, with orthonormal basis $\h{S}_k\in\mathbb{R}^{n\times r} $, be the span of the top $r$ left singular vectors of $R_k$. Ties are broken arbitrarily.

\end{itemize}
\item Return $\h{\SU}_K$.
\end{itemize}



\end{minipage}}\end{center}

\caption{$\alg$ for memory-limited PCA from incomplete data. \label{fig:Alg}}


\end{figure}

\subsection{Interpretation of $\alg$}
\label{sec:interp}
$\alg$ has a natural interpretation as a solver for a non-convex optimization program, which we next discuss. First, however, let us enrich our notation: Recall the measurement blocks $\{Y_k\}_{k=1}^K$ and let the (random) index set $\Omega_k\subseteq [1:n]\times [1:b_k]$ be the support of $Y_k$, for every $k$. We write that $Y_k=P_{\Omega_k}(S_k)$, where the columns of $S_k \in\mathbb{R}^{n\times b_k}$ belong to the true subspace $\SU$, namely $S_k = S \cdot Q_k^*$ for some coefficient matrix $Q_k\in\mathbb{R}^{b_k \times r}$. Here,  $P_{\Omega_k}(S_k)$  retains only the entries of $S_k$ on the index set $\Omega_k$, setting the rest to zero.
Note that $\{S_k\}_{k=1}^K$  and $\{Q_k\}_{k=1}^K$ are formed by partitioning $\{s_t\}_{t=1}^T$ and $\{q_t\}_{t=1}^T$ into $K$ blocks, respectively. Likewise, $\{\Omega_k\}$ are formed from $\{\omega_t\}$.



With this introduction, let us form
\begin{align}
Y := \l[
\begin{array}{cccc}
Y_1 & Y_2 & \cdots & Y_K
\end{array}
\r]
\in\mathbb{R}^{n\times T},
\end{align}
by concatenating all measurement blocks, supported on  $\Omega\subseteq [1:n]\times [1:T]$. 
To find the  true subspace $\SU$, one might then solve\begin{align}
\begin{cases}
\underset{\U,X}{\min}
%_{\U \in \mathbb{G}(n,r)}
 \l\|P_{\U^\perp}X  \r\|_F^2, \\
P_{\Omega}\l( X\r) = Y,
\end{cases}
\label{eq:non cvx}
\end{align}
where the minimization is over all $r$-dimensional subspaces $\U\subset\mathbb{R}^n$ and  matrices $X\in\mathbb{R}^{n\times T}$.
Also, $P_{\U^\perp}\in\mathbb{R}^{n\times n}$ is the orthogonal projection onto the orthogonal complement of subspace $\U$, and $P_\Omega(X)$ retains only the entries of $X$ on the index set $\Omega$, setting the rest to zero.    In particular, with full measurements ($\Omega=[1:n]\times [1:T]$), Program \eqref{eq:non cvx}  reduces to PCA, as it returns $X=Y$ and sets $\U$ equal the span of the top $r$  left singular vectors of $Y$. Note also that Program \eqref{eq:non cvx} is a non-convex problem because the Grassmannian $\mathbb{G}(n,r)$, the set of all $r$-dimensional subspaces in $\mathbb{R}^n$, is a non-convex set.\footnote{The Grassmannian can be embedded in $\mathbb{R}^{n\times n}$ via the map that takes $\U\in\mathbb{G}(n,r)$ to the orthogonal projection $P_{\U}$. The resulting submanifold of $\mathbb{R}^{n\times n}$ is a non-convex set.}

Instead of Program \eqref{eq:non cvx}, consider the equivalent program
\begin{align}
\begin{cases}
\min \sum_{k=1}^K  \l\|P_{\U_k}^\perp X_k  \r\|_F^2,\\
P_{\Omega_k}\l(X_k\r) = Y_k, \qquad  k\in[1:K],\\
\U_1 = \U_2 = \cdots = \U_K,
\end{cases}
\label{eq:non cvx parallel}
\end{align}
where the minimization is over all subspaces $\{\U_k\}\subset\mathbb{G}(n,r)$ and matrices $\{X_k\}$,  with $X_k\in\mathbb{R}^{n\times b_k}$ for every $k$.

Consider the following approximate solver for Program \eqref{eq:non cvx parallel}:
\begin{itemize}
\item Setting $X_1=Y_1$ in Program \eqref{eq:non cvx parallel},   we minimize
$ \min \|P_{\U_1}^\perp Y_1 \|_F^2
$ over $\U_1\in\mathbb{G}(n,r)$ and find a minimizer to be the span of  top $r$ left singular vectors of $Y_1$, namely $\h{\SU}_1$ in $\alg$.
\item For $k\in[2:K]$, repeat the following:
\begin{itemize}
\item Setting $\U_k=\h{\SU}_{k-1}$ in Program \eqref{eq:non cvx parallel}, we  solve
\begin{align}
\begin{cases}
\min \l\|P_{\h{\SU}_{k-1}^\perp}X_k \r\|_F^2, \\
P_{\Omega_k}\l(X_k\r) = Y_k,
\end{cases}
\label{eq:Rk interp}
\end{align}
over all matrices $X_k\in\mathbb{R}^{n\times b_k}$. In Appendix \ref{sec:interp verifying}, we verify that the minimizer of Program \eqref{eq:Rk interp} is $R_k$ in $\alg$.
\item Setting $X_k=R_k$ in Program \eqref{eq:non cvx parallel},   we minimize
$ \min \|P_{\U_k}^\perp R_k \|_F^2
$ over $\U_k\in\mathbb{G}(n,r)$ to find  $\h{\SU}_k$ in $\alg$.
\end{itemize}
\end{itemize}
Note that following this procedure precisely  produces $\{R_k\}$ and $\{\h{\SU}_k\}$  in $\alg$.

Another insight about the choice of $R_k$ in \eqref{eq:Rk} is as follows. At the beginning of the $k$th iteration of $\alg$ with $k\ge 2$, the available estimate of the true subspace is $\h{\SU}_{k-1}$ with orthonormal basis $\h{S}_{k-1}$. Given a new measurement vector $y\in\mathbb{R}^n$, supported on the index set $\omega\subseteq [1:n]$, $z= \h{S}_{k-1} (P_{\omega} \h{S}_{k-1})^\dagger y$ best approximates $y$ in $\h{\SU}_{k-1}$ (in the $\ell_2$ sense).
In order to agree with the measurements, we minimally adjust this to $y+P_{\omega^C}z$, where $P_{\omega^C}$ projects onto the complement of $\omega$. This indeed matches the expression for the columns of $R_k$ in $\alg$. We note that this type of ``least-change''  strategy has been successfully used in the development of quasi-Newton methods for optimization \cite[Chapter 6]{nocedal2006numerical}.

Lastly, it would be interesting to find a statistical interpretation for \eqref{eq:Rk}, perhaps in the form of a conditional expectation on the new measurement block. This might be a nontrivial task, however, as it seems to require forming a prior distribution on the Grassmannian.







\section{Performance of $\alg$}
\label{sec:theory}

This section summarizes the theoretical guarantees for $\alg$, with the  derivations deferred to Section \ref{sec:analysis} and the appendices.  To measure the performance of $\alg$, we naturally use  principal angles as our error metric. More specifically, recall that  $\SU$ and $\h{\SU}_K$  denote the true subspace and the output of $\alg$, respectively. Then, the $i$th largest singular value of $P_{\SU^{\perp}} P_{\h{\SU}_{K}}$ equals $\sin(\theta_i(\SU,\h{\SU}_K))$, where 
$$
\theta_1(\SU,\h{\SU}_K)\ge \theta_2(\SU,\h{\SU}_K) \ge \cdots \ge \theta_r(\SU,\h{\SU}_K)
$$
denote the principal angles between $\SU$ and $\h{\SU}_K$ \cite{golub2013matrix}.
Our error metric is then
\begin{equation}
d_{\mathbb{G}}\l(\SU,\h{\SU}_K \r) := \sqrt{
\frac{1}{r}
\sum_{i=1}^r \sin^2\l( \theta_i\l(\SU,\h{\SU}_K\r)
\r)}
=
\frac{\l\| P_{\SU^\perp}P_{\h{\SU}_K} \r\|_F}{\sqrt{r}}.
\label{eq:err metric}
\end{equation}


The concept of \emph{coherence} is also critical in  parsing our results, since we consider entrywise subsampling. The coherence of an $r$-dimensional subspace $\SU$ with orthonormal basis $S\in\mathbb{R}^{n\times r}$ is defined as
\begin{equation}
\eta\l( \SU\r) := % \frac{n}{r}  \l\| S\r\|^2_{2\rightarrow\infty} =
\frac{n}{r} \max \l\|S[i,:] \r\|_2^2,
\label{eq:def of coh}
\end{equation}
where $S[i,:]$ is the $i$th row of $S$. It is easy to verify that $\eta(\SU)$ is independent of the choice of orthonormal basis $S$, and that $1 \le \eta(\SU)\le \frac{n}{r}$.  It is also common to say that $\SU$ is \emph{coherent} (\emph{incoherent}) when $\eta(\SU)$ is large (small). Loosely speaking, when $\SU$ is coherent, its orthonormal basis $S$ is ``spiky.'' An example is when $\SU$ is the span of a column-subset of the identity matrix. In contrast, when $\SU$ is incoherent, entries of $S$ tend to be ``diffuse.'' Not surprisingly, identifying a coherent subspace from subsampled data  may require many more measurements  \cite{balzano2015local,mitliagkas2014streaming,chen2015incoherence}.

Lastly, we will often use $\lesssim$ and $\gtrsim$ to suppress universal constants, as well as the standard big-O  and big-Omega notations, $O(\cdot)$ and $\Omega(\cdot)$, respectively.  With this setup, our main result is as follows.

\begin{thm}\label{thm:main result} \textbf{\emph{[Performance of $\alg$]}}
Consider an $r$-dimensional subspace $\SU$ with orthonormal basis $S\in\mathbb{R}^{n\times r}$.  For an  integer $T$, let the coefficient vectors $\{q_t\}_{t=1}^T\subset \mathbb{R}^r$ be independent copies of a random vector $q\in\mathbb{R}^r$.
For every $t\in[1:T]$,  we observe each coordinate of  $s_t = S q_t \in \SU$ independently and with a probability of $p$,  and we collect the measurements in $y_t\in\mathbb{R}^n$, setting unobserved entries to zero.

For an integer $K$, fix block sizes $\{b_k\}_{k=1}^K$ such that $\sum_{k=1}^K b_k = T$ and $b_k\ge r$ for all $k$.
Fix also  $\alpha\gtrsim 1, \nu\ge1$ and $\{\eta_k\}_{k=1}^K$ such that  $1\le  \eta_k\le b_k/(\alpha^2 \nu^2 r)$ for all $k$.
By partitioning $\{q_t\}_t$, form $\{Q_k\}_{k=1}^K$ with $Q_k\in\mathbb{R}^{b_k\times r}$ for each $k$. Let $\nu(Q_k)$ be the condition number of $Q_k$, namely the ratio of its largest and smallest singular values. Set  also $\Q_k=\mbox{span}(Q_k)$ with coherence $\eta(\Q_k)$ (see (\ref{eq:def of coh})).

Then, with $\h{\SU}_K$  being the output of $\alg$, it holds that
\begin{equation}
d_{\mathbb{G}}\l(\SU,\h{\SU}_K \r)
   \lesssim   2^{-K}  \cdot
\alpha \cdot  \nu \cdot \sqrt{ \l( 1\vee  \frac{n}{b_1} \r) \frac{\l(\eta_1 \vee \eta\l(\SU \r) \r)
r \log(n\vee \max_k b_k) }{pn}},
\label{eq:geo decay}
\end{equation}
except with a probability of
\begin{equation}
 K\cdot O\l(e^{-\alpha}+(\min_k b_k)^{-\alpha} \r)+ \sum_{k=1}^K   \Pr\l[\nu\l( Q_{k}\r) >\nu \r] +  \Pr\l[  \eta\l(\Q_{k} \r) > \eta_{k} \r]  ,
 \label{eq:failure prob thm}
\end{equation}
and provided that
\begin{equation}
p\gtrsim \frac{\alpha^6 \nu^6}{\alpha^2\nu^2-\Cl{lwr bnd p}}  \l(1 \vee  \frac{n}{b_1}\r) \frac{\l(\eta_1\vee \eta(\SU) \r)r^2 \log^4(n\vee \max_k b_k)}{n},\qquad
p \ge 1- \frac{\Cr{lwr bnd p}}{\alpha^2 \nu^2},
\label{eq:reqs on n p}
\end{equation}
for an absolute constant $\Cr{lwr bnd p}>0$. Above, $a\vee b=\max(a,b)$ for  $a,b\in\mathbb{R}$.
\end{thm}

A few remarks are in order, starting with a simplified statement of Theorem \ref{thm:main result}. 
\begin{rem}\textbf{\emph{[Discussion]}} \label{rem:essence}\emph{Let us first discuss  the essence of Theorem \ref{thm:main result}, before making things more concrete in the next remarks.  Suppose that $\eta(\SU)=1$, namely the true subspace is incoherent. Suppose also that we partially measure ``generic'' vectors  drawn from $\SU$ and choose the following for the parameters in Theorem \ref{thm:main result}: $\alpha=O(1)$,  $\nu=O(1)$, $b_1=\Omega(n)$, $b_2=b_3=\cdots =b_K=\Omega(r\log r)$, and  $\eta_k=O(\log b_k)$ for all $k$.\footnote{Strictly speaking, $b_k$ must satisfy $r\lesssim b_k/\log b_k$ for $k\ge 2$, in light of Theorem \ref{thm:main result} and Remark \ref{rem:coeff}, specifically \eqref{eq:cnd on nu and eta}. This requirement can also be expressed by means of the \emph{Lambert W-function} \cite[\textsection 4.13]{olver2010nist}. However, $\log\log r$ is often a small constant, which justifies the simpler statement $b_k=\Omega(r\log r)$ in Remark \ref{rem:essence}. }
 In this case, Theorem \ref{thm:main result} roughly states that the estimation error of $\alg$ obeys
$$
d_{\mathbb{G}}\l(\SU,\h{\SU}_K\r) \lesssim 2^{-K} \log n \sqrt{\frac{r}{pn}},
$$
with high probability, and provided that
$$
p\gtrsim \frac{r^2\log^5 n}{n},
\qquad p\ge \Cl{display},
$$
for some constant $\Cr{display}\in(0,1)$.  In other words, the estimation error of $\alg$ decays exponentially fast provided we observe $\max(r^2\log^5n,\Cr{display}n)$ entries of every incoming vector. In the following remarks, we will discuss different aspects of this result in more details.}
\end{rem}


\begin{rem}\textbf{\emph{[Global linear convergence]}}
\emph{According to Theorem \ref{thm:main result}, $\alg$ converges globally and linearly to the true subspace, with high probability and if the sampling probability $p$ is sufficiently large. In this sense,
Theorem \ref{thm:main result} provides the first comprehensive analysis of a class of  ``least-change'' algorithms that includes $\alg$ and $\grouse$ \cite{balzano2010online}, as discussed in more detail in Section  \ref{sec:related work}.}
\end{rem}
\begin{rem} \textbf{\emph{[Sampling probability]}} \emph{A notable shortcoming of Theorem \ref{thm:main result} is the requirement in (\ref{eq:reqs on n p}) that $p\ge 1- \Cr{lwr bnd p}/\alpha^2\nu^2$.
 A perhaps less pressing issue is the quadratic (rather than linear) relation in (\ref{eq:reqs on n p}) between the ambient dimension $n$ and the typically small rank $r$. }

\emph{
 As the strong numerical performance of $\alg$ in Section \ref{sec:simulations} suggests, these drawbacks of Theorem \ref{thm:main result} appear to be artifacts of our proof technique that have resisted all efforts to remove them. See Remark \ref{rem:proof} for the technical roots of these difficulties.}
\end{rem}

\begin{rem}\label{rem:coeff} \textbf{\emph{[Coefficients]}} \emph{ Influence of the coefficient vectors $\{q_t\}_{t=1}^T$ in Theorem \ref{thm:main result}, as captured by the coherence factors $\{\eta(\Q_k)\}_{k=1}^K$, is often mild. For example, suppose that $\{q_t\}_{t=1}^T$ are independent copies of a standard Gaussian vector, in which case $\{Q_k\}_{k=1}^K$ would become independent standard random Gaussian matrices. Then, with fixed $\alpha \ge 1$ as in Theorem \ref{thm:main result}, familiar arguments in random matrix theory yield that
\begin{equation}
\nu\l(Q_k\r) = O(1),
\qquad \eta\l(Q_k\r) \lesssim \alpha \log b_k ,
\qquad \forall k\in[1:K],
\label{eq:familiar}
\end{equation}
when $\min_k b_k \gtrsim \alpha^2 r$ and except with a probability of at most $K(e^{-\alpha r}+ (\min_k b_k)^{-\alpha})$.  For the sake of completeness, (\ref{eq:familiar}) is proved in Appendix \ref{sec:cnd coh of Gaussian}.  Therefore, with the choice of
\begin{equation}
\nu=O(1),\qquad \eta_k=O(\alpha \log b_k), \qquad \forall
 k\in[1:K],
 \label{eq:cnd on nu and eta}
\end{equation}
the failure probability in (\ref{eq:failure prob thm}) reads
\begin{equation}
K\cdot O\l(e^{-\alpha}+ (\min_k b_k)^{-\alpha} \r) +
K\cdot\l( e^{-\alpha r}+ (\min_k b_k)^{-\alpha} \r)
= K\cdot O\l(e^{-\alpha}+(\min_k b_k)^{-\alpha } \r).
\end{equation}
}\emph{We close this remark by pointing out that the dependence on $\{\eta(\Q_k)\}_k$ in Theorem \ref{thm:main result} is not an artifact of the analysis. For example, suppose that   each column of every $Q_k$ contains only two nonzero entries, both equal to one, located at indexes $2i$ and $2i+1$, with integer $i$ selected uniformly at random. Then, for all $k$, $\eta(\Q_k)=\frac{n}{2r}$ is large and, even when $p=1$, it is impossible to correctly recover the true subspace $\SU$. }
\end{rem}


\begin{rem} \label{rem:block size}\textbf{\emph{[Block sizes]}} \emph{For the sake of discussion, suppose  that the coefficient vectors are independent standard Gaussian  vectors and take $\nu$ and $\{\eta_k\}_k$  as in Remark \ref{rem:coeff}.
Both (\ref{eq:geo decay}) and (\ref{eq:reqs on n p}) strongly suggest setting the first block size in $\alg$ as $b_1=\Omega (n)$. Remark \ref{rem:coeff} suggests setting the rest of block sizes as $b_k=\Omega(r)$, $k\ge 2$. The requirement on $b_1$ is likely \emph{not} an artifact, as discussed in Remark \ref{rem:proof}; there appear to be pathological scenarios in which $b_1=\Omega(n)$ is necessary for global convergence of $\alg$. However, in all  simulations in Section \ref{sec:simulations},  we  set $b_1=b_2=\cdots=b_K=\Omega(r)$ without facing any convergence problems. }

\end{rem}

\begin{rem} \label{rem:complexity} \emph{\textbf{[Computational complexity]}} \emph{We measure the algorithmic complexity of $\alg$ by calculating the average number of floating-point operations performed on  an incoming vector. Again, consider random Gaussian coefficients and the choice of $\nu$ and $\{\eta_k\}$ in Remark \ref{rem:coeff}. Then, in light of Remark \ref{rem:block size}, let us take $b_1=\Omega(n)$ and $b_k=\Omega(r)$ for $k\ge 2$. }

\emph{Iteration $k=1$ of $\alg$ involves (approximately) finding the top $r$ left singular vectors of  $Y_1\in\mathbb{R}^{n\times b_1}$ which, using a randomized SVD, could be done after $O(rn^2)$ operations or equivalently	 $O(rn)$ operations per vector \cite[Section 1.6]{halko2011finding}.}
\emph{At the $k$th iteration  with $k\ge 2$, $\alg$ requires finding the pseudo-inverse of $P_{\omega_j}\h{S}_{k-1}\in\mathbb{R}^{n\times r}$ for each incoming vector which costs $O(r^2n)$ operations. Overall, the computational complexity of $\alg$ in iteration $k\ge 2$ is $O(r^2n)$ operations  per  vector. As further discussed in Section \ref{sec:related work}, this matches the complexity of other algorithms for memory-limited PCA from limited measurements.}
\end{rem}

\begin{rem}\label{rem:storage}\emph{\textbf{[Storage]}}
\emph{We measure the storage required by $\alg$ by calculating the number of memory elements stored by $\alg$ at any given instant.  Again, consider random Gaussian coefficients and choice of $\nu$ and $\{\eta_k\}$ in Remark \ref{rem:coeff}. Then, in light of Remark \ref{rem:block size}, let us take $b_1=\Omega(n)$ and $b_k=\Omega(r)$ for $k\ge 2$.}


\emph{The fist iteration ($k=1$) of $\alg$ requires keeping $Y_1\in\mathbb{R}^{n\times b_1}$ in the memory and, therefore, has a storage requirement of $O(pn^2)$  or simply $O(n^2)$ memory elements. However, the rest of the iterations have a far lower storage requirement of $O(rn)$ memory elements.
 Indeed, at the $k$th iteration for $k\ge 2$, $\alg$ must store the current estimate $\h{S}_{k-1}\in\mathbb{R}^{n\times r}$ and the new measurement block $Y_k\in\mathbb{R}^{n\times r}$. This translates into $O(rn)+O(prn)$ or simply $O(rn)$ memory elements.
See Section \ref{sec:related work} for comparison with other algorithms.}
\end{rem}

\begin{rem}\label{rem:proof} \textbf{\emph{[Proof strategy and room for improvement]}} \emph{To prove Theorem \ref{thm:main result}, we first establish local linear convergence. That is, we establish in Lemma \ref{lem:whp main} that $\alg$ converges to the true subspace $\SU$ (with high probability and when $p$ is large enough) \emph{if} the algorithm is initialized sufficiently close to $\SU$, namely when $\h{\SU}_1$ is close enough to $\SU$. The proof of this claim crucially involves a tight perturbation bound (Lemma \ref{lem:perturbation lemma}), as well as
 controlling the spectral norms of $P_{\SU^\perp}R_k$  and $R_k - SQ_k^*$ (see (\ref{eq:Rk})), which is difficult due to the nature of $R_k$. As a proxy for the spectral norm, we bound the Frobenius norm of these quantities  in Appendix \ref{sec:Proof-of-Lemma whp main} and five ensuing appendices. This compromise, however, leads to the quadratic (rather than linear) dependence of $p$ on $r$  in (\ref{eq:reqs on n p}), which is likely an artifact of the proof technique.  Other key technical steps to obtain Lemma \ref{lem:whp main} include avoiding statistical dependence in Section \ref{sec:numerator}, as well as
 finding a tight bound on a key random variable in Lemma \ref{lem:O norm} by directly analyzing the associated ``normal equation.'' Another artifact of this argument is the constant lower bound on $p$, which appears difficult to remove.}

\emph{The next step towards the proof of Theorem \ref{thm:main result} is establishing in Lemma \ref{lem:init} that $\h{\SU}_1$, produced by $\alg$, is sufficiently close to $\SU$. This lemma allows us to guarantee that $\h{\SU}_1$ is a good initialization that would ``activate'' the local linear convergence in the rest of $\alg$ iterations ($k\ge 2$) promised by Lemma \ref{lem:whp main}. The proof of Lemma \ref{lem:init} involves a standard application of the matrix Bernstein inequality in Appendix \ref{sec:proof of lemma init}.}

\emph{Lastly, we combine Lemmas \ref{lem:init} and \ref{lem:whp main} together to complete the proof. In particular, we find that, for $\h{\SU}_1$ to be a good initialization, taking $b_1=\Omega(n)$ is necessary in the proof. Whether this requirement is an artifact of the proof is not clear to the authors. In particular, it has been recently established that any local minimum of the objective function in \eqref{eq:non cvx} (after small changes) is also a global minimum \cite{ge2016matrix}. However, it is not clear what the implication of this result is here since $\alg$ does not belong to the class of solvers in Theorem 2.3 in \cite{ge2016matrix}. Further investigating the requirement $b_1=\Omega(n)$ is needed.}
\end{rem}


\section{Related Work}
\label{sec:related work}

Among several algorithms that have been proposed for tracking low-dimensional structure in a data set from partially observed streaming measurements \cite{mitliagkas2014streaming,chi2013petrels,mardani2015subspace,xie2013change}, $\alg$ might be most closely related to $\grouse$ \cite{balzano2010online,balzano2013grouse}. $\grouse$ performs memory-limited PCA from incomplete data using stochastic gradient projection on the Grassmannian, updating  its estimate of the true subspace with each new measurement vector.

Both $\grouse$ and $\alg$ were designed based on the principal of least change, discussed in Section~\ref{sec:interp}. In fact, when $\grouse$ is sufficiently close to the true subspace and with a specific choice of its step length, both algorithms have nearly identical updates (see \cite[Equation 1.9]{balzano2015local}).



Local convergence of $\grouse$, in expectation, was recently established  in \cite{balzano2015local}. More specifically, \cite{balzano2015local} stipulates that, if the current estimate $\widehat{\SU}_k$ is sufficiently close to the true subspace $\SU$, then $\widehat{\SU}_{k+1}$ will be even closer to $\SU$ and this leads to a linear local convergence in expectation.


Our own update was inspired by that of $\grouse$ when we found zero-filled updates were unreliable. However, $\grouse$ was derived as a purely streaming algorithm, and it therefore is not designed to leverage common low-rank structure that may be revealed when a block of vectors is processed at once. Therefore for each block, $\alg$ achieves a more significant reduction in error than is possible with $\grouse$. %
 Furthermore, our guarantees hold with high probability, rather than in expectation,
 an improvement made possible by  the  machinery deployed in this work. Most importantly, $\alg$ addresses the lack of global convergence in $\grouse$ with missing data. Indeed, the first measurement block in $\alg$ provides a good initialization (see Lemma \ref{lem:init}) which, combined with Lemma \ref{lem:whp main}, guarantees global convergence of $\alg$ to the true subspace (see Theorem \ref{thm:main result}).

Lastly, both $\alg$ and $\grouse$ have a computational complexity of $O(rn)$ floating-point operations per incoming vector (see Remark \ref{rem:complexity}). Also, $\alg$ for $k\ge 2$ and $\grouse$ both require $O(rn)$ memory elements of storage. In theory, $\alg$ requires a large first block $b_1=O(n$) and consequently $O(n^2)$ memory elements of storage (see Remark \ref{rem:storage}). In all simulations in Section \ref{sec:simulations}, however, we set $b_1=O(r)$, requiring a storage of $O(rn)$ memory elements.
 With full measurements ($p=1$), a close relative of both $\alg$ and $\grouse$ is incremental SVD  \cite{bunch1978updating,balsubramani2013fast,oja1985stochastic}.






We must also discuss the algorithms in \cite{mitliagkas2014streaming,chi2013petrels,mardani2015subspace,xie2013change}.
The algorithm in \cite{mitliagkas2014streaming}, in a sense, extends the classic power method to  handle missing data, in part by improving the main result of \cite{lounici2014high}. With high probability, this algorithm converges globally and linearly to the true subspace and, most notably, succeeds for arbitrarily small sampling probability $p$, if the scope of the algorithm $T$ is large enough. Additionally, this algorithm too has a  computational complexity of   $O(rn)$ operations per vector and a storage requirement of $O(rn)$ memory elements.

$\operatorname{PETRELS}$, introduced in \cite{chi2013petrels},  operates on one column at a time (rather than blocks) and convergence is known only to a stationary point of (14), with global convergence an apparently open question. Designed for online \emph{matrix completion}, the algorithm in \cite{mardani2015subspace} also operates on one column at a time (rather than blocks) and asymptotic onvergence to the true subspace is established (see Propositions 2 and 3). This framework is also extended to tensors. $\operatorname{MOUSSE}$ in  \cite{xie2013change} tracks a union of subspaces (manifold) rather than just one;  $\alg $ would function more like an ingredient of this algorithm. Asymptotic consistency of $\operatorname{MOUSSE}$ is also established there.
We point out that the theoretical guarantees for $\alg$ surpass those of other algorithms; in Section \ref{sec:theory}, we showed nonasymptotic global convergence of $\alg$.


In the next section, we compare the performance of three of these algorithms in practice and find that $\alg$ displays state-of-the-art performance.


As discussed in Section \ref{sec:problem statement}, $\alg$ is also useful in Big Data applications where data could only be processed in small amounts. In particular, $\alg$ is closely related to {matrix completion} algorithms. For example, consider a  rank-$r$ matrix $M\in\mathbb{R}^{n\times T}$ with column span of $\SU\in\mathbb{G}(n,r)$.
Recovering $M$ from its entries on a (typically small) index set $\Omega\subset[1:n]\times [1:T]$ is known as the matrix completion problem \cite{candes2009exact,eftekhari2016mc}. If $T$ is large, then one might complete $M\in\mathbb{R}^{n\times T}$ by partitioning it into blocks of size $n\times b$, processing each block separately and in sequence, and cycling through the blocks more than once, if desired. Note that, given $\SU=\mbox{span}(M)$, completing $M$ reduces to a least-squares problem. Therefore, each iteration of $\alg$ might be adjusted to complete the corresponding block of $M$ based on the current estimate of $\SU$, cycling through the measurement blocks in sequence until convergence. Note, however, that our analysis  in Theorem \ref{thm:main result} requires statistical independence between consecutive blocks and does not apply in this form if measurement blocks are cycled more than once. We  also refer the interested reader to the discussions in on this topic  \cite{balzano2010online,mitliagkas2014streaming}.






\section{Simulations}
\label{sec:simulations}



This section consists of two parts: first, we empirically study the dependence of $\alg$ on various parameters, and second we compare $\alg$ with $\grouse$ and the algorithm in \cite{mitliagkas2014streaming}.
 In all simulations, we consider an $r$-dimensional subspace $\SU\subset\mathbb{R}^n$ and a sequence of generic vectors $\{s_t\}_{t=1}^T\subset \SU$. Every entry of these vectors is observed with probability $p\in(0,1]$ and collected in measurement vectors $\{y_t\}_{t=1}^T\subset\mathbb{R}^n$. Our objective is to estimate $\SU$ from $\{y_t\}$, as described in Section \ref{sec:problem statement}.

\paragraph{Sampling probability} We first set $n=100$, $r=5$, and let $\SU$ be a generic $r$-dimensional subspace, namely the span of an $n\times r$ standard random Gaussian matrix. For various values of probability $p$, we run $\alg$ with  block size $b=b_1=b_2=\cdots=b_K=2r=10$ and scope of $T=500r=2500$,  recording the average (over $50$ trials) estimation error $d_{\mathbb{G}}(\SU,\h{\SU}_K)$ (see \eqref{eq:err metric}). The average error versus probability is plotted in Figure \ref{fig:variable prob}.






\paragraph{Subspace dimension}
With the same setting as the previous paragraph,  we now set $p=3r/n=0.15$ and vary the subspace dimension $r$, block size $b=2r$, and scope $T=500r$. The average error versus subspace dimension is plotted in Figure \ref{fig:variable rank}.


\paragraph{Ambient dimension}
This time, we set $r=5$, $p=3r/n$, $b=2r$, $T=500r$, and vary the ambient dimension $n$. In other words, we vary $n$ while keeping the number of measurements per vector fixed at $\sim pn = 3r$.
The average error versus ambient dimension is plotted in Figure \ref{fig:variable ambient}. Observe that the performance of $\alg$ steadily degrades as $n$ increases. This is in agreement with Lemma \ref{lem:whp main} which, roughly speaking, sets the convergence rate at $\sqrt{1-p}$; as $n$ increases, $p=3r/n$ and the convergence rate both drop. See also Theorem 2.14 in \cite{balzano2015local}.

\paragraph{Block size}
Next, we set $n=100$, $r=5$, $p=3r/n$, $T=500r$, and vary the block size $b=b_1=b_2=\cdots=b_K$.
The average error versus block size in both cases is depicted in Figure \ref{fig:variable block}. As discussed in Remark \ref{rem:block size}, a block size of $b=\Omega(r)$ is necessary for the success of $\alg$, which explains the poor performance of $\alg$ for very small values of $b$. However, as $b$ increases, the number of blocks $ K=T/b$ reduces because the scope $T$ is  held fixed. As the estimation error of $\alg$ scales like $2^{-K}$ in Theorem \ref{thm:main result}, the performance suffers in Figure \ref{fig:variable block}. It appears that the choice of $b=Cr$ in $\alg$ for a relatively small  $C$ (like $C=2$) guarantees the best performance.

\paragraph{Coherence}
Lastly, we set $n=300$, $r=10$, $p=3r/n$, $b=2r$, and $T=500r$. We then test the performance of $\alg$  as the coherence of $\SU$ varies (see \eqref{eq:def of coh}). To that end, let $\SU\subset\mathbb{R}^n$ be a generic subspace with orthonormal basis $S\in\mathbb{R}^{n\times r}$. Then, the average coherence of $\SU$ over $50$ trials was $3.3334\ll n/r$ and  the average estimation error of $\alg$  was $2.795\cdot 10^{-5}$. On the other hand, let $D\in\mathbb{R}^{n\times n}$ be a diagonal matrix with entries $D[i,i]=i^{-1}$ and consider $S'=DS$. Unlike $\SU$, $\SU':=\mbox{span}(S')$ is typically a coherent subspace since the energy of $S'$ is mostly concentrated along its first few rows. This time, the average coherence of $\SU'$ over $50$ trials was $19.1773\approx n/r$ and the average estimation error of $\alg$ was $0.4286$.

\begin{figure}
\begin{center}


\subfloat[\label{fig:variable prob}]{\protect\includegraphics[width=0.49\textwidth]{variable_prob.pdf}
}
\subfloat[\label{fig:variable rank}]{\protect\includegraphics[width=0.49\textwidth]{variable_rank}
}

\subfloat[\label{fig:variable ambient}]{\protect\includegraphics[width=0.49\textwidth]{variable_ambient.pdf}
}
\subfloat[\label{fig:variable block}]{\protect\includegraphics[width=0.49\textwidth]{variable_block_fixed_scope.pdf}
}

\caption{Performance of $\alg$ as (a) sampling probability $p$, (b) subspace dimension $r$, (c) ambient dimension $n$, (d) block size $b$ vary. $\widehat{\SU}_K$ is the output of $\alg$ and $d_G(\SU,\widehat{\SU}_K)$ is its distance to the true subspace $\SU$, which generated the input of $\alg$. 
See Section \ref{sec:simulations} for details.}
\end{center}
\end{figure}




\paragraph{Comparisons}

Next, we empirically compare $\alg$ with $\grouse$ and the modified power method in \cite{mitliagkas2014streaming}. We set $n=100$, $r=5$, $T=5\cdot 10^3$, and take $\SU\subset\mathbb{R}^n$ to be a generic $r$-dimensional subspace. For $p=3r/n=0.15$, we compare the three algorithms in Figure \ref{fig:comparison}, which shows the average over $100$ trials of the estimation error of algorithms (with respect to the metric $d_{\mathbb{G}}$)  as time progresses.  For $\alg$, we used the block size of $b=2r$. Having tried to get the best performance from $\grouse$, the (diminishing) step size is set to  $ 100/t$.  For \cite{mitliagkas2014streaming}, we set the block size as $b=2n$ for the best outcome.

We also compare, for various values of the sampling probability $p$, the final estimation error of the three algorithms in Figure \ref{fig:comparison variable probability}. In both tests, $\alg$ is comparable to $\grouse$, and substantially improves over the power method. 
Note also that $\alg$, unlike many similar algorithms, has provable global convergence and, because it operates on measurement blocks, we suspect that it will show more robustness against noise. We leave that investigation to a future work.


\begin{figure}
\begin{center}
\subfloat[\label{fig:comparison}]{\protect\includegraphics[width=0.48\textwidth]{Comparison}
}


\subfloat[\label{fig:comparison variable probability}]{\protect\includegraphics[width=0.48\textwidth]{Comparison_variable_probability.pdf}
}

\caption{(a) Estimation error versus time for $\alg$, $\grouse$, and the modified power method in \cite{mitliagkas2014streaming} with a prescribed set of parameters. (b) Estimation error of these algorithms as the sampling probability $p$ varies. See Section \ref{sec:simulations} for details.}
\end{center}
\end{figure}




\section{Theory}
\label{sec:analysis}
In this section, we  prove Theorem \ref{thm:main result} in two steps. First, in Section \ref{sec:Initialization}, we  establish that the first iteration of $\alg$ provides a good initial estimate of the true subspace $\SU$. Then, in Section \ref{sec:Refinement}, we  analyse how $\alg$ iteratively refines this initial estimate.


A short word on notation is in order first. We will frequently use MATLAB's matrix notation so that, for example, $A[i,j]$ is the $[i,j]$th entry of $A$, and the row-vector $A[i,:]$ corresponds to the $i$th row of $A$. By $\{\epsilon_i\}_i \overset{\operatorname{ind.}}{\sim} \mbox{Bernoulli}(p)$, we mean that $\{\epsilon_i\}_i$ are independent Bernoulli random variables taking one with probability of $p$ and zero otherwise. Throughout, $E_{i,j}$ stands for the $[i,j]$th canonical matrix so that $E_{i,j}[i,j]=1$ is its only nonzero entry. The size of $E_{i,j}$ may be inferred from the context. As usual, $\|\cdot\|$ and $\|\cdot\|_F$ stand for the spectral and Frobenius norms. In addition, $\|A\|_\infty$ and $\|A\|_{2\rightarrow\infty}$ return the largest entry of a matrix $A$ (in magnitude) and the largest $\ell_2$ norm of the rows of $A$, respectively.

For purely aesthetic reasons, we will occasionally use the notation $a \vee b = \max(a,b)$ and $a \wedge b = \min(a,b)$.

\subsection{Initialization}
\label{sec:Initialization}


In the first iteration ($k=1$), $\alg$ simply sets $\h{\SU}_1$ to be the span of the top $r$ (left) singular vectors of  the first incoming measurement block $Y_1$. The following lemma, roughly speaking,  bounds the largest principal angle between $\SU$ and $\h{\SU}_1$ as
$$
\theta_1\l(\SU,\h{\SU}_1\r) \lesssim  \sqrt{\frac{r}{pn}},
$$
and is proved in Appendix \ref{sec:proof of lemma init} with the aid of standard large deviation bounds.
\begin{lem}
For $\alpha\ge 1$, $\nu\ge 1$, and $1 \le \eta_{1} \le \frac{b_1}{r}$, it holds that
\label{lem:init}
\begin{align}
\frac{\l\| P_{\SU^\perp}P_{\h{\SU}_1} \r\|_F}{\sqrt{r}}
 \lesssim \alpha \cdot  \nu \sqrt{
  \l( 1\vee \frac{n}{b_1} \r) \frac{\l( \eta_{1} \vee \eta\l(\SU \r)\r) r\log(n\vee b_1)}{pn}},
\label{eq:init eq}
\end{align}
except with a probability of at most $e^{-\alpha}+ \Pr[ \nu(Q_1) > \nu ]+\Pr[ \eta(\Q_1)> \eta_{1} ]$.
\end{lem}



\subsection{Refinement}
\label{sec:Refinement}


At iteration $k\in[2:K]$, $\alg$ uses the current estimate $\h{\SU}_{k-1}$ and the new measurement block $Y_k$ to produce a new estimate $\h{\SU}_k$ (of the true subspace $\SU$).  The main challenge here is to compare the  principal angles $\theta_1(\SU,\h{\SU}_k)$ and $\theta_1(\SU,\h{\SU}_{k-1})$. The following result, proved in Appendix
\ref{sec:Proof-of-Lemma whp main},  roughly states that
$$
\theta_1\l(\SU,\h{\SU}_{k}\r) \lesssim  \sqrt{1-p} \cdot \theta_1\l(\SU,\h{\SU}_{k-1}\r),
$$
provided that $\theta_1(\SU,\h{\SU}_{k-1})\lesssim \sqrt{p}$.
\begin{lem}
\label{lem:whp main}
Fix $k\in[2:K]$, $\alpha\ge1$,  $\nu\ge 1$, and $1\le \eta_{k}\le \frac{b}{r}$.
Then, except with a probability of  at most
$
O(e^{-\alpha}+b_k^{1-\alpha}) + \Pr[ \nu(Q_k) > \nu  ] +\Pr[\eta(\Q_{k}) > \eta_{k}],
$
it holds that
\begin{equation}
\frac{\left\Vert P_{\SU^{\perp}}P_{\h{\SU}_{k}}\right\Vert _{F}}{\sqrt{r}}
\lesssim
\l( \alpha\cdot \nu \left(\sqrt{\frac{\eta_{k}}{b_k}}+\sqrt{1-p} \right)
+\frac{1}{4} \r)
\frac{
\left\Vert P_{\SU^{\perp}}P_{\h{\SU}_{k-1}}\right\Vert
 _{F}}{\sqrt{r}},
 \label{eq:geo fast}
\end{equation}
provided that  $p\gtrsim\alpha^{2}\log^{3}\left(n\vee b_k\right)\cdot\eta(\widehat{\SU}_{k-1})r/n$, and
\begin{equation}
\left\Vert P_{\SU^{\perp}}P_{\h{\SU}_{k-1}}\right\Vert _{F}\lesssim  \frac{\sqrt{p}}{ \alpha \cdot \nu \log^{\frac{3}{2}} (n\vee b_k)}
.\label{eq:activate}
\end{equation}


\end{lem}




\subsection{Completing the Proof of Theorem \ref{thm:main result}}

We now combine Lemmas \ref{lem:init} and \ref{lem:whp main} to conclude the analysis of $\alg$. Fix $\alpha\ge 1$, $\nu\ge 1$, and $\{\eta_{k}\}_{k=1}^K$ with  $1\le \eta_{k}\le b_k/r$ for every $k$.
Conveniently, set $b_{\max} := \max_k b_k$. Recall  that  $\alg$ produces the sequence of estimates $\{\h{\SU}_k\}_{k=1}^K$.  In light of Lemma \ref{lem:init}, for $\h{\SU}_1$ to be a good initialization that would ``activate'' \eqref{eq:activate} in Lemma \ref{lem:whp main} with $k=2$, it suffices to have that
\begin{equation}
 \alpha \cdot  \nu \sqrt{
  \l( 1\vee \frac{n}{b_1} \r) \frac{ \l(\eta_1 \vee \eta\l(\SU \r)\r) r\log(n\vee b_1)}{pn}}
\lesssim
\frac{1}{ \alpha \cdot \nu \log^{\frac{3}{2 }} (n\vee b_2)}
\sqrt{\frac{p}{r}}.
\label{eq:1st cnd on m pre}
\end{equation}
That is, under \eqref{eq:1st cnd on m pre}, both \eqref{eq:init eq} and \eqref{eq:activate} (with $k=2$) hold except with a probability of at most $e^{-\alpha}+ \Pr[ \nu(Q_1) > \nu ]+\Pr[ \eta(\Q_1)> \eta_1 ]$. In other words, $\h{\SU}_1$ is a good initial estimate for the true subspace and, from this point onward, the estimation error reduces exponentially fast, as dictated by \eqref{eq:geo fast}. To quantify this exponential convergence, next suppose that
\begin{equation}
\label{eq:assump}
\alpha\cdot \nu \left(\sqrt{\frac{\eta_k}{b_k}}+\sqrt{1-p} \right)
\lesssim 1,
\qquad \forall k\in [1:K-1],
\end{equation}
so that, in particular, the factor in front of $\|P_{\SU^\perp}P_{\h{\SU}_{k-1}}\|_F/\sqrt{r}$ in \eqref{eq:geo fast} does not exceed $1/2$.  In particular, as long as $\eta_k\lesssim b_k/(\alpha\cdot \nu)^2$,  \eqref{eq:assump} simplifies to
\begin{equation}
p \ge 1- \frac{\Cr{lwr bnd p}}{\alpha^2 \nu^2},
\label{eq:lwr bnd simple}
\end{equation}
for an absolute constant $\Cr{lwr bnd p}>0$.  In light of \eqref{eq:lwr bnd simple}, \eqref{eq:1st cnd on m pre} holds when
\begin{equation}
 \alpha \cdot  \nu \sqrt{
  \l( 1\vee \frac{n}{b_1} \r) \frac{ \l(\eta_1 \vee \eta\l(\SU \r)\r) r\log(n\vee b_1)}{ \l(1-\frac{\Cr{lwr bnd p}}{\alpha^2\nu^2} \r) n}}
\lesssim
\frac{1}{ \alpha \cdot \nu \log^{\frac{3}{2 }} (n\vee b_2)}
\sqrt{\frac{p}{r}},
\end{equation}
or
\begin{equation}
p\gtrsim \frac{\alpha^6 \nu^6}{\alpha^2\nu^2-\Cr{lwr bnd p}} \l(1\vee \frac{n}{b_1}\r) \frac{\l(\eta_1\vee \eta(\SU) \r)r^2 \log^4(n\vee b_{\max})}{n}.
\label{eq:1st cnd on m}
\end{equation}
Assume, without loss of generality, that $\|P_{\SU^\perp}P_{\h{\SU}_k}\|_F\ne 0$ for every $k\in[1:K-1]$. Otherwise, by Lemma \ref{lem:whp main}, $\alg$ converges to $\SU$ early,  with high probability.  Then, for every $k\in[2:K]$, Lemmas \ref{lem:init} and \ref{lem:whp main} imply that
\begin{align}
& \frac{\l\| P_{\SU^\perp} P_{\h{\SU}_{k}} \r\|_F}{\sqrt{r}}
\nonumber\\
& =  \prod_{k'=2}^k \frac{\l\| P_{\SU^\perp} P_{\h{\SU}_{k'}} \r\|_F}{\l\| P_{\SU^\perp} P_{\h{\SU}_{k'-1}} \r\|_F}
\cdot  \frac{\l\| P_{\SU^\perp} P_{\h{\SU}_{1}} \r\|_F}{\sqrt{r}} \nonumber\\
& \lesssim
\prod_{k'=2}^k
\l(
\alpha\cdot \nu \left(\sqrt{\frac{\eta_k}{b}}+\sqrt{1-p} \right)+ \frac{1}{4}
  \r)
  \cdot
  \alpha \cdot  \nu \sqrt{ \l( 1\vee \frac{n}{b_1} \r) \frac{\l(\eta_1 \vee \eta\l(\SU \r) \r)
r \log(n\vee b_{\max}) }{pn}}
\qquad \mbox{(Lemmas \ref{lem:init} and \ref{lem:whp main})}
 \nonumber\\
& \lesssim   2^{-k+1}  \cdot
\alpha \cdot  \nu \sqrt{ \l( 1\vee \frac{n}{b_1} \r) \frac{\l(\eta_1 \vee \eta\l(\SU \r) \r)
r \log(n\vee b_{\max}) }{pn}},
\qquad \mbox{(see \eqref{eq:assump})}
\label{eq:iterative bnd}
\end{align}
provided that
\begin{equation}
p\gtrsim \alpha^{2}  \max_{k\in[1:K-1]}\eta(\widehat{\SU}_{k})  \cdot \frac{ r \log^{3}\left(n\vee b_{\max}\right)}{n},
\label{eq:another req on p}
\end{equation}
 and except with a probability of at most
\begin{align}
 K\cdot O\l(e^{-\alpha}+(\min_{k'} b_{k'})^{1-\alpha} \r)+ \sum_{k'=1}^K   \Pr\l[\nu\l( Q_{k'}\r) >\nu \r] +  \Pr\l[  \eta\l(\Q_{k'} \r) > \eta_{k'} \r].
%\nonumber\\
%& =   k\cdot O\l(e^{-\alpha}+b^{1-\alpha} \r)+ k \l( \Pr\l[\nu\l( Q\r) >\nu \r] +  \Pr\l[  \eta\l(\Q \r) > \eta_b \r] \r).
\label{eq:final failure prob}
\end{align}
In words, \eqref{eq:iterative bnd} states that the estimation error is halved in every iteration, with high probability and provided that \eqref{eq:lwr bnd simple},  \eqref{eq:1st cnd on m},  and \eqref{eq:another req on p} hold.

In particular, to better interpret \eqref{eq:another req on p}, we must replace the coherence $\eta(\h{\SU}_k)$ therein with a simpler quantity, perhaps  $\eta(\SU)$. We can do so thanks to Lemma \ref{lem:coherence perturbation} in the Toolbox (Appendix  \ref{sec:Toolbox}) which, roughly speaking, states that a pair of subspaces $\A$ and $\B$ with a  small principal angle  have similar coherences, namely $ \theta_1(\A,\B) \approx 0 \Longrightarrow \eta(\A)\approx \eta(\B)$. More concretely, using \eqref{eq:iterative bnd} and after invoking Lemma \ref{lem:coherence perturbation},  we find the following for every $k\in[1:K]$:
\begin{align*}
\sqrt{\eta(\h{\SU}_{k} )} & \le \sqrt{\eta\l(\SU \r)}+ \l\|P_{\SU^\perp}P_{\SU_k} \r\|  \sqrt{\frac{n}{r}}
\qquad \mbox{(see Lemma \ref{lem:coherence perturbation})}
\nonumber\\
%& = \sqrt{\eta\l(\SU\r)} + \l\|P_{\SU^\perp}P_{\SU_k} \r\| \sqrt{\frac{n}{r}}
%\qquad \mbox{(see \eqref{eq:basic id about principal angle 1})}
%\nonumber\\
& \le \sqrt{\eta\l(\SU\r)} + \l\|P_{\SU^\perp}P_{\SU_k} \r\|_F \sqrt{\frac{n}{r}} \nonumber\\
& \lesssim \sqrt{\eta\l( \SU \r)} + \sqrt{r} \l(
2^{-k+1}  \cdot
\alpha \cdot  \nu \sqrt{  \l( 1\vee \frac{n}{b_1} \r) \frac{\l(\eta_1 \vee \eta\l(\SU \r) \r)
r \log(n\vee b_{\max})}{pn}} \r)
\sqrt{\frac{n}{r}}
\qquad \mbox{(see \eqref{eq:iterative bnd})}
\nonumber\\
& = \sqrt{\eta\l(\SU\r)} + 2^{-k+1} \cdot \alpha \cdot \nu \sqrt{ \l( 1 \vee \frac{n}{b_1}\r)  \frac{\l( \eta_1\vee \eta\l(\SU\r)\r)  r
\log(n\vee b_{\max})}{p}} \nonumber\\
& \le \sqrt{\eta\l(\SU\r)} +  \alpha \cdot \nu \sqrt{ \l( 1 \vee \frac{n}{b_1}\r)  \frac{\l( \eta_1 \vee \eta\l(\SU\r)\r) r
\log(n\vee b_{\max})}{p}} \nonumber\\
& \lesssim  \alpha \cdot \nu \sqrt{ \l( 1 \vee \frac{n}{b_1}\r)  \frac{\l(\eta_1 \vee \eta\l(\SU\r) \r) r
\log(n\vee b_{\max})}{p}} \nonumber\\
& \le \alpha \cdot \nu \sqrt{ \l( 1 \vee \frac{n}{b_1}\r)  \frac{\l(\eta_1 \vee \eta\l(\SU\r) \r) r
\log(n\vee b_{\max})}{1-\Cr{lwr bnd p}/\alpha^2 \nu^2}}.
\qquad \mbox{(see \eqref{eq:lwr bnd simple})}
\end{align*}
Consequently, \eqref{eq:another req on p} holds when
\begin{align}
\label{eq:not consequential}
p \gtrsim \frac{\alpha^6 \nu^4}{\alpha^2\nu^2-\Cr{lwr bnd p}} {\l(1\vee \frac{n}{b_1}\r) \frac{\l(\eta_{1}\vee \eta(\SU)  \r) r^2 \log^4(n\vee b_{\max})}{n}}.
\end{align}
To summarize, the exponential convergence promised in \eqref{eq:iterative bnd} is valid if \eqref{eq:lwr bnd simple} and  \eqref{eq:1st cnd on m} hold (because  \eqref{eq:1st cnd on m} implies \eqref{eq:not consequential} too).
This completes the proof of Theorem \ref{thm:main result}.

\section*{Acknowledgements}
AE would like to thank Anand Vidyashankar and Chris Williams for separately pointing out the possibility of a statistical interpretation of $\alg$, as discussed at the end of Section \ref{sec:interp}. MBW was partially supported by NSF grant CCF-1409258 and NSF CAREER grant CCF-1149225.




\bibliographystyle{unsrt}
\bibliography{ReferencesSNIPE}


\appendix





\section{Toolbox}
\label{sec:Toolbox}

This section collects a number of useful results for the reader's convenience. We begin with the following concentration inequalities that are repeatedly invoked in the rest of the appendices
\cite{gross2011recovering,tropp2012user}.
\begin{lem}
\textbf{\textup{{[}Matrix Bernstein inequality for Frobenius norm{]}
}}\textup{\emph{Let $\{Z_{i}\}\subset\mathbb{R}^{n\times b}$ be a
finite sequence of zero-mean independent random matrices, and set
\[
\beta:=\max_{i}\left\Vert Z_{i}\right\Vert _{F},
\]
\[
\sigma^{2}:=\sum_{i}\mathbb{E}\left\Vert Z_{i}\right\Vert _{F}^{2}.
\]
Then, for $\alpha\ge1$ and except with a probability of at most $e^{-\alpha}$,
it holds that
\[
\left\Vert \sum_{i}Z_{i}\right\Vert _{F}\lesssim\alpha \cdot \max(\beta,\sigma).
\]
}}
\label{lem:Bernie for fro}
\end{lem}

\begin{lem}
\textbf{\emph{{[}Matrix Berstein inequality for spectral norm{]} }}\textup{\emph{Let
$\{Z_{i}\}\subset\mathbb{R}^{n\times b}$ be a finite sequence of
zero-mean independent random matrices, and set
\[
\beta:=\max_{i}\|Z_{i}\|,
\]
\[
\sigma^{2}:=\left\Vert \sum_{i}\mathbb{E}\left[Z_{i}^{*}Z_{i}\right]\right\Vert \vee\left\Vert \sum_{i}\mathbb{E}\left[Z_{i}Z_{i}^{*}\right]\right\Vert .
\]
Then, for $\alpha\ge1$ and except with a probability of at most $e^{-\alpha}$,
it holds that
\[
\left\Vert \sum_{i}Z_{i}\right\Vert \lesssim\alpha\cdot \max\left(\log\left(n\vee b\right)\cdot \beta,\sqrt{\log\left(n\vee b\right)}\cdot\sigma\right).
\]
}}
\label{lem:Bernie for spec}
\end{lem}



Consider $r$-dimensional subspaces $\A$ and $\B$.    A review of the following well-known identities is perhaps helpful:
\begin{equation}
\sin\l(\theta_1\l(\A,\B\r)\r) = \l\| P_{\A^\perp} P_{\B} \r\| = \l\| P_{\A}-P_{\B} \r\|,
\label{eq:basic id about principal angle 1}
\end{equation}
\begin{equation}
\sqrt{\sum_{i=1}^r \sin^2\l(\theta_i\l(\A,\B\r)\r)}
 =   \l\| P_{\A^\perp} P_{B} \r\|_F = \frac{1}{\sqrt{2}}\l\| P_{\A}-P_{\B} \r\|_F.
 \label{eq:basic id about principal angle 2}
\end{equation}
The following lemma is a simple variation of the standard perturbation
bounds \cite{wedin1972perturbation}.
\begin{lem}
\textbf{\emph{{[}Perturbation bound{]} }}\textup{\label{lem:perturbation lemma}
}\textup{\emph{Fix a rank-$r$ matrix $A$ and let  $\operatorname{A}=\mbox{span}(A)$.
For matrix $B$, let $B_{r}$ be a rank-$r$ truncation of $B$ obtained via SVD, and set $\operatorname{B}_{r}=\mbox{span}(B_{r})$.
Then, it holds that
\[
\left\Vert P_{\operatorname{A}}-P_{\operatorname{B}_{r}}\right\Vert =\left\Vert P_{\operatorname{A}^{\perp}}P_{\operatorname{B}_{r}}\right\Vert \le\frac{\left\Vert P_{\operatorname{A}^{\perp}}B\right\Vert }{\sigma_{r}\left(A\right)-\left\Vert B-A\right\Vert }\le\frac{\left\Vert B-A\right\Vert }{\sigma_{r}\left(A\right)-\left\Vert B-A\right\Vert },
\]
}}\textup{
\[
\frac{1}{\sqrt{2}}\left\Vert P_{\operatorname{A}}-P_{\operatorname{B}_{r}}\right\Vert _{F}=\left\Vert P_{\operatorname{A}^{\perp}}P_{\operatorname{B}_{r}}\right\Vert _{F}\le\frac{\left\Vert P_{\operatorname{A}^{\perp}}B\right\Vert _{F}}{\sigma_{r}\left(A\right)-\left\Vert B-A\right\Vert }\le\frac{\left\Vert B-A\right\Vert _{F}}{\sigma_{r}\left(A\right)-\left\Vert B-A\right\Vert },
\]
where $\sigma_{r}(A)$ is the smallest nonzero singular value of
$A$. }\end{lem}
\begin{proof}
Let $B_{r^{+}}:=B-B_{r}$ denote the residual and note that
\begin{align*}
\left\Vert P_{\A^{\perp}}P_{\B_{r}}\right\Vert  & =\left\Vert P_{\A^{\perp}}B_{r}B_{r}^{\dagger}\right\Vert \qquad\left(\B_{r}=\mbox{span}(B_{r})\right)\\
 & =\left\Vert P_{\A^{\perp}}\left(B-B_{r^{+}}\right)B_{r}^{\dagger}\right\Vert \qquad\left(B=B_{r}+B_{r^{+}}\right)\\
 & =\left\Vert P_{\A^{\perp}}BB_{r}^{\dagger}\right\Vert \qquad\left(\mbox{span}\left(B_{r^{+}}^{*}\right)\perp\mbox{span}\left(B_{r}^{*}\right)=\mbox{span}\left(B_{r}^{\dagger}\right)\right)\\
 & \le\left\Vert P_{\A^{\perp}}B\right\Vert \cdot\left\Vert B_{r}^{\dagger}\right\Vert \\
 & =\frac{\left\Vert P_{\A^{\perp}}B\right\Vert }{\sigma_{r}\left(B_{r}\right)}\\
 & \le\frac{\left\Vert P_{\A^{\perp}}B\right\Vert }{\sigma_{r}\left(A\right)-\left\Vert B-A\right\Vert }.\qquad\mbox{(Weyl's inequality)}
\end{align*}
The proof is identical for the claim with the Frobenius norm and is omitted here.
\end{proof}


\begin{lem}\textbf{\emph{[Coherence under perturbation]}}
\label{lem:coherence perturbation}
Let $\A,\operatorname{B}$ be  two $r$-dimensional subspaces in $\mathbb{R}^n$, and let $\theta_1(\A,\mbox{B})$ be their (largest) principal angle. Then, their coherences are related as follows:
\begin{equation*}
\sqrt{\eta\l( \operatorname{B} \r)}
\le \sqrt{\eta\l( \A \r)} +  \l\|P_{\A}^\perp P_{\B} \r\|\sqrt{\frac{n}{r}}.
\end{equation*}
\end{lem}

\begin{proof}
Let $\theta_1\ge \theta_2\ge \cdots \ge \theta_r$ be the principal angles between the subspaces $\A$ and $\B$ so that, in particular, $\theta_1=\theta_1(\A,\mbox{B})$. It is well-known \cite{golub2013matrix} that there exist orthonormal bases $A,B\in\mathbb{R}^{n\times r}$ for the subspaces $\A$ and $\B$, respectively, such that
\begin{equation}
A^* B = \mbox{diag}\l(
\l[
\begin{array}{cccc}
\cos \theta_1 &
\cos \theta_2 &
\cdots &
\cos \theta_r
\end{array}
\r]
  \r)
=: \Gamma \in \mathbb{R}^{r\times r},
\label{eq:cosines}
\end{equation}
where $\mbox{diag}(a)$ is the diagonal matrix formed from vector $a$.
There also exists $A'\in\mathbb{R}^{n\times r}$ with orthonormal columns such that
\begin{equation}
\l(A'\r)^*B = \mbox{diag}\l(
\l[
\begin{array}{cccc}
\sin \theta_1 &
\sin \theta_2 &
\cdots &
\sin \theta_r
\end{array}
\r]
  \r)
=: \Sigma \in\mathbb{R}^{r\times r },
\qquad \l(A'\r)^* A = 0,
\label{eq:sines}
\end{equation}
and, moreover,
\begin{equation}
\mbox{span}\l(
\l[
\begin{array}{cc}
A & B
\end{array}
\r]
 \r)
 =
 \mbox{span}\l(
 \l[
\begin{array}{cc}
A & A'
\end{array}
\r]
 \r).
 \label{eq:decomp proof}
\end{equation}
With $\A'=\mbox{span}(A')$, it follows that
\begin{align}
B & = P_{\A} B + P_{\A'} B
\nonumber\\
& = AA^* B + A' (A')^* B
\nonumber\\
& = A \Gamma + A' \Sigma.
\qquad \mbox{(see \eqref{eq:cosines} and \eqref{eq:sines})}
\label{eq:princ decomposition}
\end{align}
Consequently,
\begin{align}
\sqrt{\eta\l( \mbox{B} \r)} &
= \sqrt{\frac{n}{r}}\max_i \l\| B[i,:] \r\|_2
\qquad \mbox{(see \eqref{eq:def of coh})}
\nonumber\\
& \le  \sqrt{\frac{n}{r}}
\max_i \l\|  A[i,:] \cdot \Gamma \r\|_2 +
\sqrt{\frac{n}{r}}
\max_i \l\| A'[i,:] \cdot \Sigma  \r\|_2
\qquad \mbox{(\eqref{eq:princ decomposition} and triangle inequality)}
\nonumber\\
& \le \sqrt{\frac{n}{r}}
\max_i \l\|  A[i,:] \r\|_2 \l\| \Gamma  \r\| +
\sqrt{\frac{n}{r}}
\max_i \l\| A'[i,:] \r\|_2 \l\| \Sigma  \r\| \nonumber\\
& = \sqrt{\eta \l( \A \r)} \l\| \Gamma \r\|+
\sqrt{\eta( \A' )} \l\| \Sigma \r\|
\qquad \mbox{(see \eqref{eq:def of coh})}
\nonumber\\
& \le \sqrt{\eta \l( \A \r)} \l\| \Gamma\r\| +
\sqrt{\frac{n}{r}} \l\| \Sigma \r\|
\qquad \l(  \eta(\A ') \le \frac{n}{r}\r)
\nonumber\\
& \le \sqrt{\eta \l( \A \r)}  +
\sqrt{\frac{n}{r}} \sin \theta_1
\qquad
\mbox{(see \eqref{eq:cosines} and \eqref{eq:sines})}\nonumber\\
& = \sqrt{\eta \l( \A \r)}  +
\sqrt{\frac{n}{r}}  \l\|P_{\A^\perp}P_{\B} \r\|
\qquad \mbox{(see \eqref{eq:basic id about principal angle 1})}
\end{align}
which completes the proof of Lemma \ref{lem:coherence perturbation}.
\end{proof}





\section{Supplement to Section \ref{sec:interp}}
\label{sec:interp verifying}
In this section, we verify that $R_k$ (see \eqref{eq:Rk})  is indeed the solution of
\begin{align}
\begin{cases}
\min \l\|P_{\h{\SU}_{k-1}^\perp}X_k \r\|_F^2, \\
P_{\Omega_k}\l(X_k\r) = Y_k,
\end{cases}
\label{eq:Rk interp app}
\end{align}
when $k\ge 2$. First note that Program \eqref{eq:Rk interp app} is separable and equivalent to the following $b$ programs:
\begin{equation}
R_k[:,j] =
\arg
\begin{cases}
\min \l\| P_{\SU_{k-1}^\perp} x\r\|_2^2, \\
 P_{\omega_t} \cdot  x
 = y_t,
\end{cases}
\qquad t = (k-1)b+j,
\qquad j \in [1:b].
\label{eq:Rk separated}
\end{equation}
Here, $R_k[:,j]\in\mathbb{R}^n$ is the $j$th column of $R_k$ (in MATLAB's matrix notation). To solve each of the programs in \eqref{eq:Rk separated}, we make the change of variables $x=y_t+Q_{\omega_t^C} \cdot z$. Here, $z\in\mathbb{R}^{n-m}$ and  $Q_{\omega_t^C}\in\{0,1\}^{n\times (n-m)}$ is defined naturally so that $P_{\omega_t^C}=Q_{\omega_t^C}Q_{\omega_t^C}^*$.   We now rewrite \eqref{eq:Rk separated} as
\begin{equation}
z_j := \arg \max \l\| \l(S_{k-1}^\perp\r)^* y_t
+ \l(S_{k-1}^\perp\r)^* Q_{\omega_t^C} \cdot  z
\r\|_2^2, \\
\qquad t = (k-1)b+j,
\qquad j \in [1:b].
\end{equation}
The solutions of these least-square programs are
\begin{equation}
z_j
=
-\l( \l(S_{k-1}^\perp\r)^* Q_{\omega_t^C}\r)^\dagger \l(S_{k-1}^\perp \r)^* y_t ,\qquad j\in [1:b].
\end{equation}
For fixed $j$, we simplify the expression for $z_j$ as follows:
\begin{align}
z_j & =
-\l( \l(S_{k-1}^\perp\r)^* Q_{\omega_t^C}\r)^\dagger \l(S_{k-1}^\perp \r)^* y_t \nonumber\\
& = - \l( Q_{\omega_t^C}^* P_{\SU_{k-1}^\perp}  Q_{\omega_t^C} \r)^{-1} Q_{\omega_t^C}^* P_{\SU_{k-1}^\perp} y_t
\qquad \l(\mbox{almost surely, if } m\ge r ??\r)\nonumber\\
& = - \l( I_{n-m}- Q_{\omega_t^C}^* P_{\SU_{k-1}}  Q_{\omega_t^C} \r)^{-1} Q_{\omega_t^C}^* P_{\SU_{k-1}^\perp} y_t
\nonumber\\
& = - \l( I_{n-m} + Q_{\omega^C_t}^* S_{k-1}
\l( I_r - S_{k-1}^*P_{\omega_t^C}S_{k-1} \r)^{-1} S_{k-1}^*Q_{\omega_t^C}
  \r)
 Q_{\omega_t^C}^* P_{\SU_{k-1}^\perp} y_t
 \qquad \mbox{(inversion lemma)} \nonumber\\
 & = - \l( I_{n-m} + Q_{\omega^C_t}^* S_{k-1}
\l( I_r - S_{k-1}^*P_{\omega^C_t}S_{k-1} \r)^{-1} S_{k-1}^*Q_{\omega^C_t}
  \r)
 Q_{\omega_t^C}^* P_{\SU_{k-1}^\perp} P_{\omega_t} y_t
 \qquad \l( y_t = P_{\omega_t}y_t \r) \nonumber\\
 & =  \l( I_{n-m} + Q_{\omega_t^C}^* S_{k-1}
\l( I_r - S_{k-1}^*P_{\omega^C_t}S_{k-1} \r)^{-1} S_{k-1}^*Q_{\omega_t^C}
  \r)
 Q_{\omega_t^C}^* P_{\SU_{k-1}}  y_t
 \qquad \l( Q_{\omega_t^C} Q_{\omega_t}=0 \r) \nonumber\\
 & = Q_{\omega_t^C}P_{\SU_{k-1}}y_t
 + Q_{\omega_t^C}^* S_{k-1}
\l( I_r - S_{k-1}^*P_{\omega^C_t}S_{k-1} \r)^{-1} S_{k-1}^*P_{\omega_t^C} P_{\SU_{k-1}} y_t
\nonumber\\
 & = Q_{\omega_t^C}P_{\SU_{k-1}}y_t
 + Q_{\omega_t^C}^* S_{k-1}
\l( I_r - S_{k-1}^*P_{\omega^C_t}S_{k-1} \r)^{-1}
\l( S_{k-1}^*P_{\omega_t^C}S_{k-1}-I_r\r) S_{k-1}^*y_t
\nonumber\\
& \qquad + Q_{\omega_t^C}^* S_{k-1} \l( I_r- S_{k-1}^* P_{\omega_t^C} S_{k-1} \r)^{-1} S_{k-1}^* y_t \nonumber\\
& = Q_{\omega_t^C}^* S_{k-1} \l( I_r- S_{k-1}^* P_{\omega_t^C} S_{k-1} \r)^{-1} S_{k-1}^* y_t \nonumber\\
& = Q_{\omega_t^C}^* S_{k-1} \l( S_{k-1}^* P_{\omega_t} S_{k-1} \r)^{-1} S_{k-1}^* y_t \nonumber\\
& = Q_{\omega_t^C}^* S_{k-1} \l( P_{\omega_t}S_{k-1}\r)^{\dagger}  y_t,
\qquad \l( y_t = P_{\omega_t}y_t\r)
\end{align}
which means that
\begin{equation}
y_t + Q_{\omega_t^C} z_j =
y_t + P_{\omega_t^C} S_{k-1}\l(P_{\omega_t}S_{k-1}   \r)^{\dagger} y_t,
\end{equation}
is the solution of the $j$th program in  \eqref{eq:Rk separated} which matches the $j$th column of $R_k$ defined in \eqref{eq:Rk}.







\section{Proof of Lemma \ref{lem:init}}
\label{sec:proof of lemma init}



Consider the  coefficient vectors $\{q_t\}_{t=1}^{b_1}\subset\mathbb{R}^r$ and $\{s_t=S q_t\}_{t=1}^{b_1}\subset\mathbb{R}^n$, constructed in Section \ref{sec:problem statement}. By concatenating $\{q_t\}$ and $\{s_t\}$, we form $Q_1\in\mathbb{R}^{b_1\times r}$ and $S_1=SQ_1^*\in\mathbb{R}^{n\times b_1}$, respectively. Each $s_t$ is observed on the random index set $\omega_t\subseteq[1:n]$ or, equivalently, $S_1$ is observed on the random index set $\Omega_1\subseteq [1:n]\times [1:b_1]$.
 We write this measurement process as $Y_1=P_{\Omega_1}(S_1)$, where $P_{\Omega_1}(\cdot)$ projects onto index set $\Omega_1$.


Let us fix $Q_1$ for now. Also let $Y_{1,r}\in \mathbb{R}^{n\times b_1}$ be a rank-$r$ truncation of $Y_1$, obtained via SVD. $\alg$ then  sets $\h{\SU}_1=\mbox{span}(Y_{1,r})$. Our objective here is to control $\|P_{\SU^\perp} P_{\h{\SU}_1}\|_F$. Since
\begin{equation}
\l\|P_{\SU^\perp} P_{\h{\SU}_1}\r\|_F \le \sqrt{r} \l\|P_{\SU^\perp} P_{\h{\SU}_1}\r\|,
\qquad \l( \h{\SU}_1\in\mathbb{G}(n,r) \r)
\label{eq:rank-r adv}
\end{equation}
it suffices to bound the spectral norm (instead of the Frobenius norm). Conditioned on $Q_1$, it is easy to  verify that $\mathbb{E}[Y_1]=\mathbb{E}[P_{\Omega_1}(S_1)]= p\cdot S_1$, suggesting that we might consider $Y_1$ as a perturbed copy of $p\cdot S_1$, and perhaps consider  $\h{\SU}_1=\mbox{span}(Y_{1,r})$ as a perturbation of $\SU=\mbox{span}(p\cdot S_1)$. Indeed, Lemma \ref{lem:perturbation lemma} in Toolbox (Appendix \ref{sec:Toolbox}) dictates that
\begin{align}
\l\| P_{\SU^{\perp}} P_{\h{\SU}_1} \r\|
& \le \frac{\l\| Y_1-p  S_1 \r\|}{p\cdot \sigma_r\l(   S_1 \r)-\l\| Y_1-p S_1 \r\|} \nonumber\\
& = \frac{\l\| Y_1-p S_1 \r\|}{p \cdot \sigma_r\l( Q_1 \r)-\l\| Y_1-p S_1 \r\|}
\qquad \l( S_1 = S Q_1^*,\,\, S^*S =I_r \r)
\nonumber\\
& \le \frac{2}{p}\cdot  \frac{\l\| Y_1-p S_1 \r\|}{\sigma_r\l( Q_1\r)}.
\qquad \l(   \mbox{if } \l\| Y_1-p S_{Q_1}\r\|\le  \frac{p}{2} \cdot \sigma_r(Q_1)\r)
\label{eq:perturb init}
\end{align}
It remains to bound the norm in the last line above.
To that end,
we study the concentration of $Y_1$ about its expectation by writing that
\begin{align}
Y_1 - p  S_1
&= P_{\Omega_1} \l(S_1\r) -p S_{1}
\nonumber\\
& = \sum_{i,j} \l(\epsilon_{i,j}-p\r)  S_1[i,j]\cdot E_{i,j}
\nonumber\\
& =: \sum_{i,j } Z_{i,j},
\label{eq:Bernie prep initialization}
\end{align}
where  $\{\epsilon_{i,j}\}\overset{\operatorname{ind.}}{\sim} \mbox{Bernoulli}(p)$, and $E_{i,j}\in\mathbb{R}^{n\times b_1} $ is the $[i,j]$th canonical matrix.  Additionally, $\{Z_{i,j}\}$ are independent zero-mean random matrices. In order to appeal to the matrix  Bernstein inequality (Lemma \ref{lem:Bernie for spec}), we compute the $\beta$ and $\sigma$ parameters below, starting with $\beta$:
\begin{align}
\l\| Z_{i,j} \r\| & = \l\| \l(\epsilon_{i,j}-p\r) S_1[i,j]\cdot E_{i,j} \r\| \nonumber\\
& = \l| \l(\epsilon_{i,j}-p\r) S_1[i,j] \r|\qquad  \l( \l\| E_{i,j} \r\|=1\r)\nonumber\\
& \le \l|S_1[i,j] \r| \qquad \l( \epsilon_{i,j}\in\{0,1\}  \r) \nonumber\\
& \le \l\| S_{1}\r\|_\infty
\nonumber\\
& \le \l\| S\r\|_{2\rightarrow\infty} \l\|Q_1 \r\|_{2\rightarrow\infty}
\qquad \l(S_{1} = SQ_1^*,\,\, \|AB^*\|_\infty\le \|A\|_{2\rightarrow\infty} \|B\|_{2\rightarrow\infty} \r)\nonumber\\
& \le \sqrt{\frac{\eta\l(\SU \r)r}{n}}\cdot
\sqrt{\frac{\eta\l( \Q_1 \r)r}{b_1}} \cdot \|Q_1\|
\qquad \l( \mbox{see \eqref{eq:def of coh}}
\r) \nonumber\\
& =: \beta.
\label{eq:beta in init}
\end{align}
Above, $\|A\|_\infty$ and $\|A \|_{2\rightarrow\infty}$, respectively, return the largest entry of $A$ (in magnitude) and the largest $\ell_2$ norm of the rows of matrix $A$. As for $\sigma$, we write that
\begin{align}
\l\| \mathbb{E} \l[ \sum_{i,j} Z_{i,j} Z_{i,j}^* \r] \r\|&
 = \l\| \sum_{i,j} \E \l[ \l(\epsilon_{i,j}-p \r)^2 \r] S_1[i,j]^2 \cdot E_{i,i} \r\|
\nonumber\\
& = \l\| \sum_{i,j} p(1-p) S_1[i,j]^2 \cdot E_{i,i} \r\|
\qquad \l( \epsilon_{i,j} \sim \mbox{Bernoulli}(p) \r) \nonumber\\
 & \le p \l\|\sum_{i,j} S_1[i,j]^2 \cdot E_{i,j}\r\| \nonumber\\
 & = p \l\| \sum_i \l\|S_1[i,:]\r\|^2_2 \cdot E_{i,i} \r\| \nonumber\\
 & = p \max_i \l\| S_1[i,:] \r\|_2^2 \nonumber\\
 & = p \l\| S_1 \r\|_{2\rightarrow\infty}^2 \nonumber\\
 & \le p \l\|S \r\|_{2\rightarrow\infty}^2 \cdot \l\| Q_1\r\|^2
 \qquad \l( S_1 = SQ_1^*, \,\,  \|AB\|_{2\rightarrow\infty} \le \|A\|_{2\rightarrow\infty}\|B\|\r)
\nonumber\\
& = p \cdot \frac{\eta\l( \SU\r)r}{n} \cdot \|Q_1\|^2.
\qquad \l(\mbox{see \eqref{eq:def of coh}}\r)
\label{eq:sigma leg 1 spec norm}
\end{align}
In a similar fashion, we find that
\begin{align}
\l\| \mathbb{E} \l[ \sum_{i,j} Z_{i,j}^* Z_{i,j} \r] \r\|
& \le  p  \l\| \sum_j \l\| S_1[:,j]\r\|_2^2 E_{j,j} \r\| \nonumber\\
& = p\l\| S_1^* \r\|_{2\rightarrow\infty}^2 \nonumber\\
& \le p\cdot \|S\|^2 \cdot \l\|Q_1\r\|_{2\rightarrow\infty}^2
\qquad \l(S_{1}=SQ_1^*,\,\, \l\| AB\r\|_{2\rightarrow\infty} \le \|A\|_{2\rightarrow\infty} \|B\|  \r)
\nonumber\\
& \le p \cdot \frac{\eta\l( \Q_1\r)r}{b_1} \cdot \|Q_1\|^2, \qquad \l(\|S\|=1,\,\,\mbox{see \eqref{eq:def of coh}} \r)
\label{eq:sigma leg 2 spect norm}
\end{align}
and, eventually,
\begin{align}
\sigma^2 & = \l\| \mathbb{E}\l[ \sum_{i,j} Z_{i,j}^*Z_{i,j}\r]\r\|
\vee \l\| \mathbb{E}\l[ \sum_{i,j} Z_{i,j} Z_{i,j}^*\r]\r\|
\nonumber\\
& \le \frac{pr}{n}  \l( 1\vee \frac{n}{b_1} \r)
\l( \eta\l( \SU\r)\vee \eta\l(\Q_1 \r)\r)\|Q_1\|^2.
\qquad \l( \mbox{see \eqref{eq:sigma leg 1 spec norm} and \eqref{eq:sigma leg 2 spect norm}}\r)
% \nonumber\\
%& \le  \frac{pr}{n}  \l( 1\vee \frac{n}{b} \r)
% \eta\l( \SU\r)\eta\l(\Q_1 \r) \|Q_1\|^2
% \qquad \l( \eta\l(\SU\r) \ge 1\r)
\label{eq:sigma in initialization}
\end{align}
Lastly,
\begin{align}
& \max\l( \log(n\vee b_1) \cdot \beta , \sqrt{\log(n\vee b_1)} \cdot \sigma \r)
\nonumber\\
&
\lesssim \max\l(  \log(n\vee b_1)\cdot  \frac{r}{n} ,
\sqrt{\log(n\vee b_1)} \cdot
 \sqrt{\frac{{pr}}{n}}
 \r)
\sqrt{ 1\vee \frac{n}{b_1}}
\cdot
\sqrt{\eta\l(\SU\r) \vee  \eta\l( \Q_1\r)}\cdot  \|Q_1\|
\qquad \mbox{(see \eqref{eq:beta in init} and \eqref{eq:sigma in initialization})}
\nonumber\\
& \le \sqrt{\log(n\vee b_1)} \cdot \sqrt{\frac{{pr}}{n}}
\sqrt{ 1\vee \frac{n}{b_1}} \cdot
\sqrt{\eta\l(\SU\r)\vee  \eta\l( \Q_1\r)}\cdot  \|Q_1\|.
\qquad \l(\mbox{if } p\ge \frac{\log\l( n\vee b_1\r)r}{n} \r)
\label{eq:Bernie init another prep }
\end{align}
The Bernstein inequality now dictates that
\begin{align}
\l\| Y_1 - p S_{1} \r\| & = \l\| \sum_{i,j} Z_{i,j}\r\|
\qquad \mbox{(see \eqref{eq:Bernie prep initialization})}
 \nonumber\\
& \lesssim \alpha  \max\l(\log(n\vee b_1)\cdot \beta,
\sqrt{\log(n\vee b_1} \cdot \sigma  \r)
\qquad \mbox{(see Lemma \ref{lem:Bernie for spec})}
\nonumber\\
& \lesssim
\alpha \sqrt{\log(n\vee b_1)} \cdot \sqrt{\frac{{rp}}{n}} \sqrt{1 \vee \frac{n}{b_1} }
\cdot
\sqrt{\eta\l(\SU \r) \vee \eta\l(\Q_1 \r)} \cdot \|Q_1\|,
\qquad \mbox{(see \eqref{eq:Bernie init another prep })}
\label{eq:norm bound in init}
\end{align}
except with a probability of at most $e^{-\alpha}$.
 In particular, suppose that
\begin{equation}
 p \gtrsim \alpha^2      \nu(Q_1)^2
  \l( 1\vee \frac{n}{b_1} \r)\frac{\l(\eta\l( \SU \r) \vee \eta\l( \Q_1\r) \r)r \log(n\vee b_1) }{n},
\end{equation}
so that \eqref{eq:perturb init} holds. Then, by  substituting \eqref{eq:norm bound in init} back into \eqref{eq:perturb init} and then applying \eqref{eq:rank-r adv}, we find that
\begin{align}
\frac{\l\| P_{\SU^\perp}P_{\h{\SU}_1} \r\|_F}{\sqrt{r}} &
\le  \l\| P_{\SU^\perp} P_{\h{\SU}_1} \r\| \qquad \mbox{(see \eqref{eq:rank-r adv})}
\nonumber\\
& \le \frac{2}{p} \cdot
\frac{\l\| Y_1-pS_{1} \r\|}{\sigma_{r}(Q_1)}
\qquad \mbox{(see \eqref{eq:perturb init})} \nonumber\\
& \lesssim \alpha  \sqrt{\log(n\vee b_1) \cdot \frac{r}{pn} \l( 1\vee \frac{n}{b_1} \r)  \l( {\eta\l(\SU \r)\vee \eta\l( \Q_1\r) }\r) }\cdot \nu\l( Q_1\r)
\qquad \mbox{(see \eqref{eq:norm bound in init})}\nonumber\\
& =: \delta_1\l( \nu(Q_1),\eta(Q_1)\r),
\label{eq:bnd on frob iter 1 pre final}
\end{align}
except with a probability of at most $e^{-\alpha}$, and for fixed $Q_1$. In order to remove the conditioning on $Q_1$, fix $\nu\ge 1$, $1\le \eta_{1} \le \frac{b_1}{r}$, and recall the following inequality for events $\mathcal{A}$ and $\mathcal{B}$:
\begin{align}
\Pr\l[ \mathcal{A} \r]
& = \Pr\l[ \mathcal{A} | \mathcal{B}\r] \cdot \Pr\l[ \mathcal{B}\r] + \Pr\l[ \mathcal{A} | \mathcal{B}^C\r] \cdot \Pr\l[ \mathcal{B}^C\r]\nonumber\\
& \le \Pr\l[ \mathcal{A} | \mathcal{B}\r]  +\Pr\l[ \mathcal{B}^C\r].
\label{eq:useful event ineq}
\end{align}
Set $\Q_1=\mbox{span}(Q_1)$ and let $\mathcal{E}$ be the event where both $\nu(Q_1) \le \nu$ and $\eta(\Q_1)\le \eta_{1}$.
Thanks to the inequality above, we find that
\begin{align}
& \Pr\l[ \frac{\l\|P_{\SU^{\perp}} P_{\h{\SU}_1}\r\|_F}{\sqrt{r}}
\gtrsim
 \delta_1\l( \nu, \eta_{1}\r)
\r]
 \nonumber\\
&
\le \Pr\l[ \frac{\l\|P_{\SU^{\perp}} P_{\h{\SU}_1}\r\|_F}{\sqrt{r}}
\gtrsim
\delta\l( \nu,\eta_{1}\r) \mid
\mathcal{E}
\r]
+ \Pr\l[ \mathcal{E}^C \r]
\qquad \mbox{(see \eqref{eq:useful event ineq})}
\nonumber\\
& \le e^{-\alpha} + \Pr\l[ \nu(Q_1) > \nu \r]+ \Pr\l[ \eta(Q_1)> \eta_{1} \r],
\qquad \mbox{(see \eqref{eq:bnd on frob iter 1 pre final})}
\end{align}
which completes the proof of Lemma \ref{lem:init}.

\section{Proof of Lemma \ref{lem:whp main} \label{sec:Proof-of-Lemma whp main}}

\subsection{Setup and Notation}
Let us first simplify our notation for the rest of this section. At the start of the iteration $k\in[2:K]$, we have available an estimate $\SU_o=\h{\SU}_{k-1}$ of the true subspace $\SU$ (with $o$ standing for $old$). In a sense, $\SU_o$ is our prior knowledge at this point about $\SU$. To unburden the notation, we write $Q=Q_k$ and so forth, thereby dropping the subscripts. By construction,  the columns of the coefficient matrix $Q\in\mathbb{R}^{b\times r}$ are independent copies of the random vector $q$.
Setting   $S_Q:=S\cdot Q^*\in\mathbb{R}^{n\times b}$, we  observe each entry of $S_Q$ independently with a probability of $p$, collect the measurements in $Y\in\mathbb{R}^{n\times b}$, and let $\Omega\subseteq [1:n]\times [1:b]$ denote the corresponding random index set (over which we observe $S_Q$). We write this as $Y=P_{\Omega}(S_Q)$, where $P_\Omega(\cdot)$ projects onto the index set $\Omega$.


Given the new measurement block $Y$, we update our estimate of the subspace $\SU$  (from the old $\SU_{o}=\h{\SU}_{k-1}$ to the new $\SU_n=\h{\SU}_k$) as follows: Calculate the random matrix
\begin{equation}
\mathbb{R}^{n\times b}\ni R = Y+ P_{\Omega^{C}}\left(O\left(Y\right)\right),\label{eq:main}
\end{equation}
where $P_{\Omega^C}(\cdot)$ projects onto the complement of index set $\Omega$, and
\begin{equation}
O(Y)=\left[\cdots S_{o}\cdot\left(P_{\omega_{j}}\cdot S_{o}\right)^{\dagger}\cdot Y[:,j]\cdots\right]\in\mathbb{R}^{n\times b}.
\label{eq:def of O recall}
\end{equation}
As usual, $S_o\in\mathbb{R}^{n\times b}$ is an orthonormal basis for the subspace $\SU_o$.
Assume that $Q$ is fixed for now.
Let $R_{r}$ denote a rank-$r$ truncation of $R$, obtained via SVD. Then, our updated estimate is $\SU_{n}=\mbox{span}(R_{r})$.

To control the (largest) principal angle $\theta_1(\SU,\SU_n)$, our strategy is to consider $R$ as a perturbed copy of $S_{Q}=SQ^{*}$ and, in turn,  $\SU_{n}=\mbox{span}(R_{r})$ as a perturbed copy of  $\SU=\mbox{span}(S_{Q})$.
Indeed, an application of Lemma \ref{lem:perturbation lemma}
yields that
\begin{align}
\left\Vert P_{\SU^{\perp}}P_{\SU_{n}}\right\Vert _{F}& \le\frac{\left\Vert P_{\SU^{\perp}}R\right\Vert _{F}}{\sigma_r\left(S_Q\right)-\left\Vert R-S_{Q}\right\Vert }\nonumber\\
& = \frac{\left\Vert P_{\SU^{\perp}}R\right\Vert _{F}}{\sigma_r\left(Q\right)-\left\Vert R-S_{Q}\right\Vert }
.
\qquad \l(S_Q = SQ^*,\,\, S^*S=I_r  \r)
\label{eq:objective}
\end{align}
We next bound the numerator and denominator in the last line above in separate subsections.
\subsection{Bounding the Numerator of \eqref{eq:objective}}
\label{sec:numerator}
To control the numerator, we begin with some
preparation. First, recalling the definition of $O(Y)$ from \eqref{eq:def of O recall}, we observe that
\begin{align}
O\left(Y\right) & = O\l( P_\Omega\l(S_Q\r) \r)
\qquad \l(Y=P_\Omega\l( S_Q\r)  \r)
\nonumber\\
 & =O\left(P_\Omega\l( P_{\SU_{o}}S_{Q}\right) \r) +O\left(P_\Omega\l( P_{\SU_{o}^{\perp}}S_{Q} \r)\right)\qquad\mbox{(linearity)}\nonumber \\
 & =P_{\SU_{o}}S_{Q}+O\left(P_\Omega\l(P_{\SU_{o}^{\perp}}S_{Q}\r)\right).\qquad\mbox{(see (\ref{eq:def of O recall}))}\label{eq:basic decomp}
\end{align}
The above decomposition allows us to rewrite $R$ in \eqref{eq:main} as
\begin{align}
R & = Y+ P_{\Omega^{C}}\left(O\left(Y\right)\right)\qquad\mbox{(see (\ref{eq:main}))}\nonumber \\
& =  P_\Omega(S_Q)+ P_{\Omega^{C}}\left(O\left(Y\right)\right)\qquad
\l( Y=P_\Omega(S_Q) \r)\nonumber\\
 & =P_{\Omega}\left(S_{Q}\right)+P_{\Omega^{C}}\left(P_{\SU_{o}}S_{Q}\right)+\left[P_{\Omega^{C}}\circ O\circ P_\Omega\right]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)
 \qquad\left(\mbox{see (\ref{eq:basic decomp}), }f\circ g(x):=f(g(x))\right)\nonumber \\
 & =S_{Q}-P_{\Omega^{C}}\left(S_{Q}\right)+P_{\Omega^{C}}\left(P_{\SU_{o}}S_{Q}\right)+\left[P_{\Omega^{C}}\circ  O \circ P_\Omega \right]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)
 \nonumber \\
 & = S_{Q}-P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)+\left[P_{\Omega^{C}}\circ O \circ P_\Omega \right]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)
 \nonumber \\
 & =S_{Q}-P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)+\l[ O\circ P_\Omega  \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)-\left[P_{\Omega}\circ O\circ P_\Omega\right]\left(P_{\SU_{o}^{\perp}}S_{Q}\right),\label{eq:exp for R}
\end{align}
and, consequently,
\[
P_{\SU^{\perp}}R= -
P_{\SU^\perp} \cdot
P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)+
P_{\SU^\perp} \cdot
\l[ O \circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)-
P_{\SU^\perp} \cdot
\left[P_{\Omega}\circ O \circ P_\Omega \right]\left(P_{\SU_{o}^{\perp}}S_{Q}\right).
\qquad
\l( P_{\SU^{\perp}}S_{Q}=P_{\SU^{\perp}}SQ^{*}=0
\r)
\]
In particular, with an application of the triangle inequality, it immediately  follows that
\begin{align}
\left\Vert P_{\SU^{\perp}}R\right\Vert _{F} &
\le
\left\Vert P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}+\left\Vert P_{\SU^{\perp}}\cdot \l[ O \circ P_\Omega\r] \left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}+\left\Vert \left[P_{\Omega}\circ O \circ P_\Omega \right]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}.\label{eq:2nd decomp}
\end{align}
We proceed by controlling each norm on the right-hand side above through
a series of technical lemmas (proved in Appendices \ref{sec:Proof-of-Lemma t1}-\ref{sec:Proof-of-Lemma t3}).
\begin{lem}
\label{lem:t1}For $\alpha\ge1$ and except with a probability of at
most $e^{-\alpha}$, it holds that
\[
\left\Vert P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}\lesssim\alpha\max\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}},\sqrt{1-p}\right)\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}\cdot\|Q\|.
\]
\end{lem}
\begin{lem}
\label{lem:t2}For $\alpha\ge1$ and except with a probability of at
most $b^{1-\alpha}$, it holds that
\[
\left\Vert P_{\SU^{\perp}}\cdot \l[ O\circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}\lesssim {\frac{\alpha\log^{\frac{3}{2}}\left(n\vee b\right)}{\sqrt{p}}}\cdot \left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert \cdot\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}\cdot\|Q\|,
\]
provided that $p\gtrsim\alpha^{2}\log^{3}\left(n\vee b\right)\cdot\eta(\SU_{o})r/n$.
\end{lem}

\begin{lem}
\label{lem:t3} For $\alpha\ge1$ and except with a probability of at
most $e^{-\alpha}$, it holds that
\[
\left\Vert \l[ P_{\Omega}\circ O \circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}
\lesssim
\alpha\max\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}},\sqrt{1-p}\right) \left\Vert  P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}\cdot\|Q\|.
\]
\end{lem}
Applying Lemmas \ref{lem:t1}-\ref{lem:t3}  to (\ref{eq:2nd decomp}) yields that
\begin{align}
&\left\Vert P_{\SU^{\perp}}R\right\Vert _{F}\nonumber\\
& \le
\left\Vert P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}+\left\Vert P_{\SU^{\perp}}\cdot \l[ O\circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}+\left\Vert \left[P_{\Omega}\circ O \circ P_\Omega \right]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}
\qquad \mbox{(see \eqref{eq:2nd decomp})} \nonumber\\
&
\lesssim
\alpha
\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}}+\sqrt{1-p}+{\frac{\log^{\frac{3}{2}}\left(n\vee b\right)\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert}{\sqrt{p}}} \right)\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}\cdot\|Q\|,\label{eq:numerator}
\end{align}
except with a probability of at most $O(e^{-\alpha}+b^{1-\alpha})$, and provided
that $p\gtrsim\alpha^{2}\log^{3}\left(n\vee b\right)\cdot\eta(\SU_{o})r/n$.

\subsection{Controlling the Denominator of \eqref{eq:objective}}
We now find an upper bound
for $\|R-S_{Q}\|$ in \eqref{eq:objective}. We do so in Appendix \ref{sec:t4} and summarize the outcome below.
\begin{lem}
\label{lem:t4}For $\alpha\ge1$ and except with a probability of at
most $b^{1-\alpha}$, it holds that
\[
\left\Vert R-S_{Q}\right\Vert \le \left\Vert R-S_{Q}\right\Vert _{F}\lesssim \frac{\alpha \log^{\frac{3}{2}}\left(n\vee b\right)}{\sqrt{p}}\cdot  {\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}} \cdot \|Q\| ,
\]
provided that $p\gtrsim \alpha\log^{3}\left(n\vee b\right)\cdot\eta(\SU_{o})r/n $.
\end{lem}
In particular, if we assume that
\begin{equation}
\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}\lesssim
\frac{\sqrt{p}}{\alpha \cdot  \nu(Q)\log^{\frac{3}{2}}\left(n\vee b\right)},\label{eq:initial}
\end{equation}
then Lemma \ref{lem:t4} implies that
\begin{equation}
\left\Vert R-S_{Q}\right\Vert \le\frac{\sigma_{r}(Q)}{2},\label{eq:simple denom}
\end{equation}
with high probability and when $p$ is large enough, as Lemma \ref{lem:t4} stipulates.

\subsection{Completing the Proof of Lemma \ref{lem:whp main}}
For fixed $Q$,  combining (\ref{eq:numerator}) and \eqref{eq:simple denom} allows us to finally bound the expression in \eqref{eq:objective} and arrive at
\begin{align}
 & \left\Vert P_{\SU^{\perp}}P_{\SU_{n}}\right\Vert _{F}\nonumber \\
 & \le\frac{\left\Vert P_{\SU^{\perp}}R\right\Vert _{F}}{\sigma_{r}\left(Q\right)-\left\Vert R-S_{Q}\right\Vert }
\qquad \mbox{(see \eqref{eq:objective})}
 \nonumber \\
 & \lesssim \alpha\cdot \nu(Q)\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}}+\sqrt{1-p}+{\frac{ \log^{\frac{3}{2}}\left(n\vee b\right)\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert}{\sqrt{p}}} \right)\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F},
 \label{eq:dragon pre}
\end{align}
except with a probability of at most $O(e^{-\alpha}+b^{1-\alpha})$, under (\ref{eq:initial}), and provided
that $p\gtrsim\alpha^{2}\log^{3}\left(n\vee b\right)\cdot\eta(\SU_{o})r/n$.
%and provided that the factor in front of $\|P_{\SU^\perp}P_{\SU_o}\|$   does not exceed one.
If we further assume that
\begin{equation}
\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert
\le \frac{\sqrt{p}}{4\alpha\cdot \nu(Q)\log^{\frac{3}{2}}(n\vee b) }
,\label{eq:initial 2}
\end{equation}
then \eqref{eq:dragon pre} simplifies to
\begin{align}
\left\Vert P_{\SU^{\perp}}P_{\SU_{n}}\right\Vert _{F}
&
\lesssim \l( \alpha \cdot \nu(Q) \left(\sqrt{\frac{\eta\left(\Q\right)r}{b}}+\sqrt{1-p} \right)
+ \frac{1}{4} \r)
\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}\nonumber\\
%& \lesssim  \alpha\cdot \nu(Q) \log^{2}\left(n\vee b\right)\cdot\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}}+\sqrt{\frac{\eta\left(\SU_{o}\right)r}{n}}+\sqrt{1-p} \right)
%\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}
%\nonumber\\
& =: \delta_2\l( \nu(Q),\eta\l(\Q\r) \r).
\label{eq:dragon}
\end{align}
In summary, conditioned on $Q$ and except with a probability of at most $O(e^{-\alpha}+b^{1-\alpha})$, (\ref{eq:dragon}) holds if (\ref{eq:initial}),(\ref{eq:initial 2})
 are met and
 $p\gtrsim\alpha^{2}\log^{3}\left(n\vee b\right)\cdot\eta(\SU_{o})r/n$.
% the factor in front of $\|P_{\SU^\perp}P_{\SU_o}\|$  in \eqref{eq:dragon} does not exceed one.
 In order to remove the conditioning on $Q$, fix $\nu \ge 1$ and $1\le \eta\le \frac{b}{r}$. Let also $\mathcal{E}$ be the event where both $\nu(Q) \le \nu$ and $\eta(\Q)\le \eta $.
Then, we apply \eqref{eq:useful event ineq} to find that
\begin{align}
& \Pr\l[ \l\|P_{\SU^{\perp}} P_{\SU_n}\r\|_F
\gtrsim
\delta_2\l( \nu,\eta\r)
\r]
 \nonumber\\
&
\le \Pr\l[ \l\|P_{\SU^{\perp}} P_{\SU_n}\r\|_F
\gtrsim
\delta_2\l( \nu,\eta\r)\,\mid\,  \mathcal{E}\r]
+ \Pr\l[ \mathcal{E}^C \r]
\qquad \mbox{(see \eqref{eq:useful event ineq})}
\nonumber\\
& \le O\l( e^{-\alpha}+b^{1-\alpha} \r) + \Pr\l[ \mathcal{E}^C \r] \nonumber\\
& \le O\l( e^{-\alpha}+b^{1-\alpha} \r) + \Pr\l[\nu(Q)>\nu\r] + \Pr\l[\eta(\Q)>\eta\r].
\qquad \mbox{(union bound)}
\end{align}
This completes the proof
of Lemma \ref{lem:whp main}.



\section{Proof of Lemma \ref{lem:t1}\label{sec:Proof-of-Lemma t1}}

With $Q$ fixed, note that
\begin{align*}
\mathbb{E}\left[P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right] & =
\mathbb{E}\l[\sum_{i,j} \l(1-\epsilon_{i,j}\r) \cdot \l( P_{\SU_o^\perp}S_Q \r)[i,j] \cdot E_{i,j} \r]
\nonumber\\
& =
(1-p)\cdot P_{\SU_{o}^{\perp}}S_{Q},
\qquad \l(\epsilon_{i,j} \sim \mbox{Bernoulli}(p) \r)
\end{align*}
where $\{\epsilon_{i,j}\}_{i,j}\overset{\operatorname{ind.}}{\sim} \mbox{Bernoulli}(p)$, and  $E_{i,j}\in\mathbb{R}^{n\times b}$ is the $[i,j]$th canonical
matrix. We write the deviation from this expectation as
\begin{align}
P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)-(1-p)\cdot P_{\SU_{o}^{\perp}}S_{Q} & =\sum_{i,j}\l( \l(1-\epsilon_{i,j}\r)- \l(1-p \r) \r) \cdot \left(P_{\SU_{o}^{\perp}}S_{Q}\right)[i,j]\cdot E_{i,j}
\nonumber\\
 & =:\sum_{i,j}Z_{i,j},
 \label{eq:def of Zs t1}
\end{align}
where  $\{Z_{i,j}\}\subset\mathbb{R}^{n\times b}$ are independent
and zero-mean random matrices. In order to appeal to the Bernstein inequality, we first compute the $\beta$ parameter:
\begin{align}
 \left\Vert Z_{i,j}\right\Vert _{F}
 & =\left\Vert \l(p-\epsilon_{i,j} \r) \cdot \left(P_{\SU_{o}^{\perp}}S_{Q}\right)[i,j]\cdot E_{i,j}\right\Vert _{F}
\qquad \mbox{(see \eqref{eq:def of Zs t1})}
 \nonumber\\
& = \l| \l( p-\epsilon_{i,j}\r) \cdot \l(P_{\SU_o^\perp} S_Q\r)[i,j]   \r|
\qquad \l( \l\| E_{i,j}\r\|_F=1 \r) \nonumber\\
& \le \l|\l(P_{\SU_o^\perp}S_Q \r)[i,j]  \r|
\qquad \l( \epsilon_{i,j}\in\{0,1\} \r) \nonumber\\
  & \le\left\Vert P_{\SU_{o}^{\perp}}S_{Q}\right\Vert _{\infty}\nonumber\\
 & \le \l\|  P_{\SU_o^\perp} S\r\|_{2\rightarrow\infty} \cdot \l\| Q \r\|_{2\rightarrow\infty}
\qquad \l( S_Q = SQ^*,\,\, \l\|AB^*\r\|_{\infty}\le \|A\|_{2\rightarrow\infty}\cdot  \|B\|_{2\rightarrow\infty} \r)
 \nonumber\\
  & \le  \l\| P_{\SU_o^{\perp}} P_{\SU} \r\|_{2\rightarrow\infty}
 \cdot \sqrt{\frac{\eta\l( \Q \r) r}{b}} \|Q\|
\qquad \l(  \mbox{see \eqref{eq:def of coh}}\r)
 \nonumber\\
%  & \le2\left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{2\rightarrow\infty}\cdot\left\Vert Q\right\Vert \qquad\left(S_{Q}=SQ^{*},\,\|AB\|_{\infty}\le\|A\|_{2\rightarrow\infty}\cdot\|B\|\right)
% \nonumber\\
% & =: 2 \sqrt{ \frac{\eta\l(\SU,\SU_o \r)}{n}} \l\|P_{\SU_o^\perp}S\r\|_F\cdot \|Q\|
&  =:\beta.
\label{eq:b pre lemma t1}
\end{align}
Above, $\|A\|_{\infty}$ and $\|A\|_{2\rightarrow\infty}$ return
the largest entry of $A$ (in magnitude) and the largest $\ell_{2}$ norm of the
rows of $A$, respectively.  As for the $\sigma$ parameter, we write that
\begin{align}
\mathbb{E}\l[ \sum_{i,j}\left\Vert Z_{i,j}\right\Vert _{F}^{2}\r] &
= \sum_{i,j} \mathbb{E}\l[\l( p-\epsilon_{i,j} \r)^2 \r] \cdot \l|\l( P_{\SU_o^\perp} S_Q \r)[i,j] \r|^2
\qquad \mbox{(see \eqref{eq:def of Zs t1})} \nonumber\\
& = p(1-p) \cdot  \sum_{i,j} \l|\l( P_{\SU_o^\perp} S_Q \r)[i,j] \r|^2
\qquad \l( \epsilon_{i,j} \sim \mbox{Bernoulli}(p)  \r)
\nonumber\\
& = p(1-p) \cdot \l\|P_{\SU_o^\perp}S_Q \r\|_F^2 \nonumber\\
 & \le (1-p) \cdot \left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{F}^{2}\cdot\|Q\|^{2}\qquad\left(S_{Q}=SQ^{*},\quad\|AB\|_{F}\le\|A\|_{F}\cdot\|B\|\right)\nonumber\\
 & =:\sigma^{2},
 \label{eq:sigma pre lemma t1}
\end{align}
and, finally,
\begin{align*}
\max\l(\beta, \sigma\r) & \lesssim\max\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}},\sqrt{1-p}\right)\cdot\left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{F}\cdot\|Q\|.
\qquad \mbox{(see \eqref{eq:b pre lemma t1} and \eqref{eq:sigma pre lemma t1})}
\end{align*}
Using the Bernstein inequality and for
arbitrary $\alpha\ge1$, it follows that
\begin{align}
\left\Vert P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)-(1-p)\cdot P_{\SU_{o}^{\perp}}S_{Q}\right\Vert _{F} & =\left\Vert \sum_{i,j}Z_{i,j}\right\Vert _{F}
\qquad \mbox{(see \eqref{eq:def of Zs t1})}
\nonumber\\
& \lesssim \alpha \max\l(\beta,\sigma \r)
\qquad \mbox{(see Lemma \ref{lem:Bernie for fro})}
 \nonumber\\
 & \lesssim\alpha\max\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}},\sqrt{1-p}\right)\cdot\left\Vert P_{\SU_o^{\perp}}S\right\Vert _{F}\cdot\|Q\|,
 \label{eq:t1 pre final}
\end{align}
except with a probability of at most $e^{-\alpha}$. Consequently, with an application of the triangle inequality, we find that
\begin{align*}
\left\Vert P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F} &
\lesssim (1-p)\left\Vert P_{\SU_{o}^{\perp}}S_{Q}\right\Vert _{F}+\alpha\max\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}},\sqrt{1-p}\right)\left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{F}\|Q\|
\qquad \mbox{(triangle inequality and \eqref{eq:t1 pre final})}
\\
 & \le (1-p)\left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{F}\cdot\|Q\|+\alpha\max\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}},\sqrt{1-p}\right)\left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{F}\|Q\|
\qquad \l(S_Q = SQ^* \r)
 \\
 & \lesssim\alpha\max\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}},\sqrt{1-p}\right)\left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{F}\|Q\|\nonumber\\
  & =\alpha\max\left(\sqrt{\frac{\eta\left(\Q\right)r}{b}},\sqrt{1-p}\right)\left\Vert P_{\SU_{o}^{\perp}}P_{\SU}\right\Vert _{F}\|Q\|,
  \qquad \mbox{(rotational invariance of Frobenius norm)}
\end{align*}
which completes the proof of Lemma \ref{lem:t1}.


\section{Proof of Lemma \ref{lem:t2}\label{sec:Proof-of-Lemma t2}}

Throughout, we fix the coefficient matrix $Q$. We begin by noting
that
\begin{align}
\left\Vert P_{\SU^{\perp}}\cdot \l[ O \circ P_\Omega  \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F} & =\left\Vert P_{\SU^{\perp}}P_{\SU_{o}}\cdot \l[ O\circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}\qquad\mbox{(see (\ref{eq:def of O recall}))}\nonumber \\
 & \le\left\Vert P_{\SU^{\perp}}P_{\SU_{o}}\right\Vert \cdot\left\Vert \l[O \circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}\qquad\left(\|AB\|_{F}\le\|A\|\cdot\|B\|_{F}\right).\label{eq:brk .39}
\end{align}
In Appendix \ref{sec:Proof-of-Lemma O}, we bound the random norm in the last line above.
\begin{lem}
\label{lem:O norm}For $\alpha\ge1$ and except with a probability of
at most $2b^{1-\alpha}$, it holds that
\[
\left\Vert \l[ O\circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}\lesssim {\frac{\alpha\log^{\frac{1}{2}}n\log b}{\sqrt{p}}}\cdot\left\Vert P_{\SU^{\perp}}P_{\SU_{o}}\right\Vert _{F}\cdot\|Q\|,
\]
provided that $p\gtrsim\alpha{}^{2}\log n\log^{2}b\cdot\eta(\SU_{o})r/n$.
\end{lem}
We therefore conclude
that
\begin{align}
\left\Vert P_{\SU^{\perp}}\cdot \l[O \circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}
& \le\left\Vert P_{\SU^{\perp}}P_{\SU_{o}}\right\Vert \cdot\left\Vert \l[ O \circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}
\qquad \mbox{(see \eqref{eq:brk .39})}
\nonumber\\
& \lesssim {\frac{\alpha\log^{\frac{1}{2}}n\log b}{\sqrt{p}}}\left\Vert P_{\SU^{\perp}}P_{\SU_{o}}\right\Vert \cdot\left\Vert P_{\SU^{\perp}}P_{\SU_{o}}\right\Vert _{F}\cdot\|Q\|,
\qquad \mbox{(see Lemma \ref{lem:O norm})}
\end{align}
except with a probability of at most $2b^{1-\alpha}$ and provided
that $p\gtrsim\alpha{}^{2}\log n\log^{2}b\cdot\eta(\SU_{o})r/n$. This
completes the proof of Lemma \ref{lem:t2}.


\section{Proof of Lemma \ref{lem:t3} \label{sec:Proof-of-Lemma t3}}

Throughout, $Q$ is fixed. Recalling the definition of $O(\cdot)$ from \eqref{eq:def of O recall},
we write that
\begin{align}
\left\Vert \l[ P_{\Omega}\circ O \circ P_\Omega\r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}^{2} & =\sum_{j=1}^{b}\left\Vert \left(P_{\omega_{j}}S_{o}\right)\left(P_{\omega_{j}}S_{o}\right)^{\dagger}\cdot
P_{\omega_j} \l( P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\r) \right\Vert _{2}^{2}\qquad\mbox{(see (\ref{eq:def of O recall}))}\nonumber \\
& =\sum_{j=1}^{b}\left\Vert \left(P_{\omega_{j}}S_{o}\right)\left(P_{\omega_{j}}S_{o}\right)^{\dagger}\cdot
 P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*} \right\Vert _{2}^{2}\qquad
\l( \l( P_{\omega_j}S_o\r)^\dagger = \l( S_o^*P_{\omega_j}S_o\r)^{-1} S_o^* P_{\omega_j}\r)
 \nonumber \\
 & =\sum_{j=1}^{b}\left\Vert P_{\SU_{o,j}}\cdot P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right\Vert _{2}^{2}\qquad\left(\SU_{o,j}:=\mbox{span}\left(P_{\omega_{j}}S_{o}\right)\right).
 \label{eq:brk 9.99}
 \end{align}
Let also $\SU_{o,j}^C:=\mbox{span}(P_{\omega_j^C}S_o)$,  and note that $S_o = P_{\omega_j}S_o+P_{\omega_j^C}S_o$ and $(P_{\omega_j}S_o)^* (P_{\omega_j^C}S_o)=0$. Consequently, $\SU_{o,j} \perp \SU_{o,j}^C$  and $P_{\SU_o}=P_{\SU_{o,j}}+P_{\SU_{o,j}^C}$, from which it follows that
\begin{align}
\l\| P_{\SU_{o,j}} \cdot P_{\SU_o^\perp} S \cdot Q[j,:]^* \r\| & =  \l\| \l( P_{\SU_o}-P_{\SU_{o,j}^C} \r) \cdot P_{\SU_o^\perp} S \cdot Q[j,:]^* \r\| \nonumber\\
& = \l\| P_{\SU_{o,j}^C}  \cdot P_{\SU_o^\perp} S \cdot Q[j,:]^* \r\|.
\label{eq:brk 9.999}
\end{align}
Using the above identity to simplify \eqref{eq:brk 9.99} leads to
 \begin{align}
 \left\Vert \l[ P_{\Omega}\circ O \circ P_\Omega\r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}^{2} & =\sum_{j=1}^{b}\left\Vert P_{\SU_{o,j}} \cdot P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right\Vert _{2}^{2}
\qquad \mbox{(see \eqref{eq:brk 9.99})}
 \nonumber \\
 & = \sum_{j=1}^{b}\left\Vert P_{\SU_{o,j}^C} \cdot P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right\Vert _{2}^{2}
\qquad \mbox{(see \eqref{eq:brk 9.999})}
 \nonumber\\
 & = \sum_{j=1}^{b}\left\Vert P_{\SU_{o,j}^C} P_{\omega_j^C}\cdot P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right\Vert _{2}^{2}
\qquad
\l( P_{\SU_{o,j}^C} = P_{\omega_j^C}S_o \l( S_o^* P_{\omega_j^C} S_o\r)^{-1} S_o^* P_{\omega_j^C}  \r)
 \nonumber\\
 & \le \sum_{j=1}^{b}\left\Vert  P_{\omega_j^C}\cdot P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right\Vert _{2}^{2}
 \nonumber\\
 & = \l\| P_{\Omega^C}\l(P_{\SU_o^\perp}SQ^* \r) \r\|_F^2
 \nonumber\\
 & = \l\| P_{\Omega^C}\l(P_{\SU_o^\perp}S_Q \r) \r\|_F^2
 \qquad \l(S_Q = SQ^* \r)
 \nonumber\\
 &
\lesssim \alpha^2\max\left({\frac{\eta\left(\Q\right)r}{b}},{1-p}\right) \left\Vert  P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}^2\cdot\|Q\|^2,
\qquad \mbox{(see Lemma \ref{lem:t1})}
\label{eq:brk 10}
\end{align}
except with a probability of at most $e^{-\alpha}$.  This completes the proof of Lemma \ref{lem:t3}.

\section{Proof of Lemma \ref{lem:t4}\label{sec:t4}}

After recalling (\ref{eq:exp for R}) and for $\alpha\ge1$, we observe that
\begin{align*}
& \left\Vert R-S_{Q}\right\Vert _{F} \nonumber\\
& \le \left\Vert -P_{\Omega^{C}}\left(P_{\SU_{o}^{\perp}}S_{Q}\right)+\left[P_{\Omega^{C}}\circ O \circ P_\Omega \right]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}\qquad\mbox{(see (\ref{eq:exp for R}))}\\
 & \le \left\Vert P_{\SU_{o}^{\perp}}S_{Q}\right\Vert _{F}+ \left\Vert \l[O \circ P_\Omega \r] \left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}
\qquad \l(\mbox{triangle inequality, } P_{\Omega^C}(\cdot) \mbox{ is a non-expansive operator} \r)
 \\
 & \le \left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{F}\cdot\|Q\|+ \left\Vert \l[ O\circ P_\Omega \r]\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}\quad\left(S_{Q}=SQ^{*},\quad\|AB\|_{F}\le\|A\|_{F}\cdot\|B\|\right)\\
 & \lesssim \left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{F}\cdot\|Q\|+ \frac{\alpha\log^{\frac{1}{2}}n\log b}{\sqrt{p}}\cdot \left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}\cdot\|Q\|\qquad\mbox{(see Lemma \ref{lem:O norm})}\\
 & \lesssim \frac{\alpha\log^{\frac{3}{2}}\left(n\vee b\right)}{\sqrt{p}}\cdot\left\Vert P_{\SU^{\perp}}P_{\SU_o}\right\Vert _{F}\cdot\|Q\|,
\end{align*}
except with a probability of at most $b^{1-\alpha}$, and provided
that $p\gtrsim\alpha\log^{3}\left(n\vee b\right)\cdot\eta\left(\SU_{o}\right)r/n$.
This completes the proof of Lemma \ref{lem:t4}.




\section{Proof of Lemma \ref{lem:O norm}\label{sec:Proof-of-Lemma O}}

Using the definition of $O(\cdot)$ in \eqref{eq:def of O recall}, we write that
\begin{align}
\left\Vert \l[O \circ P_\Omega\r] \left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}^{2} &
 =\sum_{j=1}^{b}\left\Vert S_{o}\left(P_{\omega_{j}}S_{o}\right)^{\dagger}P_{\omega_j}\left(P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right)\right\Vert _{2}^{2}
 \qquad\left(\mbox{(\ref{eq:def of O recall}) and }S_{Q}=SQ^{*}\right)\nonumber \\
 & =\sum_{j=1}^{b}\left\Vert S_{o}\left(P_{\omega_{j}}S_{o}\right)^{\dagger}\left(P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right)\right\Vert _{2}^{2}\qquad\left(
\l( P_{\omega_j}S_o\r)^{\dagger} = \l(S_o^* P_{\omega_j} S_o \r)^{-1} S_o^* P_{\omega_j}
 \right)\nonumber \\
 & =\sum_{j=1}^{b}\left\Vert \left(P_{\omega_{j}}S_{o}\right)^{\dagger}\left(P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right)\right\Vert _{2}^{2}.\qquad\left(S_{o}^{*}S_{o}=I_{r}\right)\label{eq:brk .4}
\end{align}
For fixed $j\in [1:b]$, consider the summand in the last line above:
\begin{align}
\left\Vert \left(P_{\omega_{j}}S_{o}\right)^{\dagger}\left(P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right)\right\Vert _{2} & \le\left\Vert \left(P_{\omega_{j}}S_{o}\right)^{\dagger}P_{\SU_{o}^{\perp}}\right\Vert \cdot\left\Vert P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right\Vert _{2}\nonumber \\
 & =\left\Vert \left(P_{\omega_{j}}S_{o}\right)^{\dagger}S_{o}^{\perp}\right\Vert \cdot\left\Vert P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right\Vert _{2}\qquad\left(P_{\SU_{o}^{\perp}}= S_o^{\perp} \l(S_o^{\perp}\r)^* \right) \nonumber \\
 & =:\left\Vert \widehat{Z}_{j}\right\Vert \cdot\left\Vert P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right\Vert _{2}.\label{eq:brk .5}
\end{align}
Above, as usual, $S_o^\perp$ is an orthonormal basis for the subspace $\SU_o^\perp$. We can now revisit (\ref{eq:brk .4}) and write that
\begin{align}
\left\Vert O\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F}^{2} & =\sum_{j=1}^{b}\left\Vert \left(P_{\omega_{j}}S_{o}\right)^{\dagger}\left(P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right)\right\Vert _{2}^{2}\qquad\mbox{(see (\ref{eq:brk .4}))}\nonumber \\
 & \le\max_{j}\left\Vert \widehat{Z}_{j}\right\Vert ^{2}\cdot\sum_{j=1}^{b}\left\Vert P_{\SU_{o}^{\perp}}S\cdot Q[j,:]^{*}\right\Vert _{2}^{2}\qquad\mbox{(see (\ref{eq:brk .5}))}\nonumber \\
 & =\max_{j}\left\Vert \widehat{Z}_{j}\right\Vert ^{2}\cdot\left\Vert P_{\SU_{o}^{\perp}}SQ^{*}\right\Vert _{F}^{2}\nonumber \\
 & \le\max_{j}\left\Vert \widehat{Z}_{j}\right\Vert ^{2}\cdot\left\Vert P_{\SU_{o}^{\perp}}S\right\Vert _{F}^{2}\|Q\|^{2}.\qquad\left(\|AB\|_{F}\le\|A\|_{F}\cdot\|B\|\right)\label{eq:brk .51}
\end{align}
It remains to control the maximum in the last line above. First, we focus on controlling $\|\widehat{Z}_{j}\|$ for
fixed $j\in[1:b]$. Observe
that $\widehat{Z}_{j}$ is a solution of the least-squares problem
\[
\widehat{Z}_{j}:=\arg\min_{Z\in\mathbb{R}^{n\times(n-r)}}\left\Vert S_{o}^{\perp}-\left(P_{\omega_{j}}S_{o}\right)Z\right\Vert _{F}^{2},
\]
and, therefore, satisfies the \emph{normal equation}
\[
\left(P_{\omega_{j}}S_{o}\right)^{*}\left(\left(P_{\omega_{j}}S_{o}\right)\widehat{Z}_{j}-S_{o}^{\perp}\right)=0,
\]
which is itself equivalent to
\begin{equation}
\left(S_{o}^{*}P_{\omega_{j}}S_{o}\right)\widehat{Z}_{j}=S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}.\qquad\left(P_{\omega_{j}}^{2}=P_{\omega_{j}}\right)\label{eq:norma}
\end{equation}
In fact, since
\[
\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}\right]=p\cdot S_{o}^{*}S_{o}^{\perp}=0,
\]
\[
\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}\right]=p\cdot I_{r},
\qquad \l( S_o^* S_o = I_r \r)
\]
we can rewrite (\ref{eq:norma}) as
\[
\left(S_{o}^{*}P_{\omega_{j}}S_{o}-\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}\right]\right)\widehat{Z}_{j}+p\cdot\widehat{Z}_{j}=S_{o}^{T}P_{\omega_{j}}S_{o}^{\perp}-\mathbb{E}\left[S_{o}^{T}P_{\omega_{j}}S_{o}^{\perp}\right].
\]
An application of the triangle inequality immediately implies that
\begin{align}
p\left\Vert \widehat{Z}_{j}\right\Vert  & \le\left\Vert S_{o}^{*}P_{\omega_{j}}S_{o}-\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}\right]\right\Vert \cdot\left\Vert \widehat{Z}_{j}\right\Vert +\left\Vert S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}-\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}\right]\right\Vert .\label{eq:key}
\end{align}
To control $\|\widehat{Z}_{j}\|$, we therefore need to derive large
devation bounds for the two remaining norms on the right-hand side
above. For the first spectral norm, we write that
\begin{equation}
\left\Vert S_{o}^{*}P_{\omega_{j}}S_{o}-\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}\right]\right\Vert =\left\Vert \sum_{i}\l(\epsilon_i-p \r)\cdot S_{o}^{*}E_{i,i}S_{o}\right\Vert =:\left\Vert \sum_{i}A_{i}\right\Vert ,
\label{eq:def of As}
\end{equation}
where $\{\epsilon_i\}_i \overset{\operatorname{ind.}}{\sim} \mbox{Bernoulli}(p)$, and  $E_{i,i}\in\mathbb{R}^{n\times n}$ is the $[i,i]$th canonical
matrix. Furthermore, $\{A_{i}\}_i\subset\mathbb{R}^{r\times r}$ are independent
and zero-mean random matrices. To apply the Bernstein inequality, we first compute $\beta$ as
\begin{align}
\left\Vert A_{i}\right\Vert  & =\left\Vert \l(\epsilon_i-p \r)\cdot S_o^*E_{i,i}S_o \right\Vert
\qquad \mbox{(see \eqref{eq:def of As})}
\nonumber\\
 & \le \left\Vert S_{o}^{*}E_{i,i}S_{o}\right\Vert  \qquad
\l( \epsilon_{i}\in \{0,1\} \r)
 \nonumber\\
 & = \left\Vert S_{o}[i,:]\right\Vert _{2}^{2}
 \nonumber\\
& \le \frac{\eta\left(\SU_{o}\right)r}{n}=:\beta.\qquad \qquad \mbox{(see \eqref{eq:def of coh})}
\label{eq:beta for Ais}
\end{align}
To compute $\sigma$, we write that
\begin{align}
\l\| \E\l[  \sum_{i} A_i^2 \r]  \r\|
&
= \l\| \sum_i \mathbb{E}\l[\l(\epsilon_i-p\r)^2 \r]  \l( S_o^* E_{i,i}S_o\r)^2 \r\|
\qquad \mbox{(see \eqref{eq:def of As})}
\nonumber\\
& = \l\|  \sum_i p(1-p) \l( S_o^* E_{i,i} S_o\r)^2  \r\|
\qquad \l(\epsilon_i \sim \mbox{Bernoulli}(p)  \r) \nonumber\\
& \le p \l\|\sum_i \l(S_o^* E_{i,i}S_o\r)^2  \r\| \nonumber\\
& = p \l\| \sum_i S_o^* E_{i,i} S_oS_o^* E_{i,i} S_o \r\| \nonumber\\
& = p \l\| \sum_i \l\|S_o[i,:] \r\|_2^2 \cdot S_o^* E_{i,i} S_o \r\|
\nonumber\\
& \le p \cdot \max_i \l\| S_o[i,:] \r\|_2^2 \cdot \l\| \sum_i S_o^* E_{i,i} S_o  \r\|  \nonumber\\
& = p \l\| \sum_i S_o^* E_{i,i}   S_o \r\|
\nonumber\\
& = p \cdot \frac{\eta\l( \SU_o\r)r}{n} \cdot \l\| \sum_i S_o^* E_{i,i} S_o\r\|
\qquad \mbox{(see \eqref{eq:def of coh})}
\nonumber\\
& = p \cdot \frac{\eta\l( \SU_o\r)r}{n}
\qquad \l( \sum_i E_{i,i}=I_n,\,\, S_o^* S_o = I_r \r)
\nonumber\\
& =: \sigma^2.
 \label{eq:sigma for Ais}
\end{align}
It also follows that
\begin{align}
\max\left(\log r\cdot \beta,\sqrt{\log r}\cdot\sigma\right)&
=\max\left(\frac{\log r\cdot\eta\left(\SU_{o}\right)r}{n},\sqrt{\frac{{\log r\cdot p\cdot \eta\left(\SU_{o}\right)r}}{n}}\right)
\qquad \mbox{(see \eqref{eq:beta for Ais} and  \eqref{eq:sigma for Ais})} \nonumber\\
&
\le\sqrt{\frac{{\log r\cdot p \cdot \eta\left(\SU_{o}\right)r}}{n}}.\qquad\left(\mbox{if }p \ge\frac{\log r \cdot \eta(\SU_{o})r}{n}\right)
\label{eq:max b sigma Ais}
\end{align}
As a result, for $\alpha\ge1$ and except with a probability of at
most $e^{-\alpha}$, it holds that
\begin{align}
\left\Vert S_{o}^{*}P_{\omega_{j}}S_{o}-\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}\right]\right\Vert & \lesssim\alpha\max\left(\log r\cdot \beta,\sqrt{\log r}\cdot\sigma\right)
\qquad \mbox{(see Lemma \ref{lem:Bernie for spec})}
\nonumber\\
& \le \alpha \sqrt{\frac{{\log r \cdot p \cdot \eta\left(\SU_{o}\right)r}}{n}}.\label{eq:large dev 1}
\end{align}
On the other hand, in order to apply the Bernstein inequality to the
second spectral norm in (\ref{eq:key}), we write that
\begin{equation}
\left\Vert S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}-\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}\right]\right\Vert  =\left\Vert \sum_{i} \l(\epsilon_i-p \r) S_{o}^{*}E_{i,i}S_{o}^{\perp}\right\Vert =:\left\Vert \sum_{i}A_{i}\right\Vert ,
\label{eq:def of other As}
\end{equation}
where $\{\epsilon_i \}_i\overset{\operatorname{ind.}}{\sim} \mbox{Bernoulli}(p)$, $E_{i,i}\in \mathbb{R}^{n\times n}$ is the $i$th  canonical  matrix, and  $\{A_{i}\}_i\subset\mathbb{R}^{r\times(n-r)}$ are zero-mean and
independent random matrices. To compute the $\beta$ parameter, we write that
\begin{align}
\left\Vert A_{i}\right\Vert  & =\left\Vert
\l(\epsilon_i-p  \r)
 S_{o}^{*}E_{i,i}S_{o}^{\perp}\right\Vert
\qquad \mbox{(see \eqref{eq:def of other As})}
  \nonumber \\
  & \le \l\|  S_o^* E_{i,i} S_o^\perp \r\| \qquad
  \l( \epsilon_i \in\{0,1\} \r) \nonumber\\
 & \le\left\Vert S_{o}^{*}E_{i,i}\right\Vert
\qquad \l( \l( S_o^\perp\r)^* S_o^\perp = I_{n-r} \r)
\nonumber\\
 & =\left\Vert S_{o}[i,:]\right\Vert _{2}\nonumber\\
 & \le\sqrt{\frac{\eta\left(\SU_{o}\right)r}{n}}=:\beta.
 \qquad \mbox{(see \eqref{eq:def of coh})}
\label{eq:b Ais 2}
\end{align}
To compute the $\sigma$ parameter, we notice that
\begin{align}
\l\| \E\l[ \sum_{i} A_i A_i^* \r] \r\| &
= \l\|   \sum_i \mathbb{E}\l[ \l( \epsilon_i-p\r)^2 \r] S_o^* E_{i,i} S_o^\perp \l( S_o^\perp \r)^* E_{i,i} S_o \r\|
\qquad \mbox{(see \eqref{eq:def of other As})}
\nonumber\\
& = \l\| \sum_i p(1-p) \cdot S_o^* E_{i,i} S_o^\perp \l( S_o^\perp \r)^* E_{i,i} S_o \r\|
\qquad \l( \epsilon_i \sim \mbox{Benoulli}(p) \r) \nonumber\\
& \le p \l\| \sum_i S_o^* E_{i,i} S_o^\perp \l( S_o^\perp \r)^* E_{i,i} S_o  \r\| \nonumber\\
& \le p
\l\| \sum_i S_o^* E_{i,i}  E_{i,i} S_o  \r\|
\qquad \l(  S_o^\perp \l( S_o^\perp \r)^* \preceq I_n  \r)
\nonumber\\
& = p \l\| \sum_i S_o^* E_{i,i}  S_o  \r\|  \nonumber\\
& = p. \qquad \l( \sum_i E_{i,i}=I_n,\,\, S_o^* S_o = I_r \r)
 \label{eq:pre sigma Ais leg 1}
\end{align}
In a similar fashion, we find that
\begin{align}
\left\Vert \mathbb{E}\l[\sum_{i}A_{i}^{*}A_{i}\r]\right\Vert
& \le p \l\|\sum_i \l(S_o^\perp \r)^* E_{i,i} S_o S_o^* E_{i,i} S_o^\perp   \r\|
\nonumber\\
&= p \l\| \sum_i \l\| S_o[i,:]\r\|_2^2 \cdot   \l(S_o^\perp \r)^* E_{i,i} S_o^\perp \r\| \nonumber\\
& \le p \cdot \max_i \l\|S_o[i,:]  \r\|_2^2 \cdot \l \| \sum_i \l(S_o^\perp \r)^* E_{i,i} S_o^\perp \r\| \nonumber\\
& = p \cdot \max_i \l\|S_o[i,:]  \r\|_2^2
\qquad \l( \sum_i E_{i,i}=I_n,\,\, \l(S_o^\perp \r)^* S_o^\perp = I_{n-r} \r) \nonumber\\
& = p \cdot \frac{\eta\l(\SU_o\r)r}{n},
\qquad \mbox{(see \eqref{eq:def of coh})}
 \label{eq:pre sigma Ais leg 2}
\end{align}
and, finally,
\begin{align}
\sigma
& =
\max\l( \l\| \E \sum_{i} A_i A_i^*\r\|,
\l\| \E \sum_{i} A_i^*A_i  \r\|
   \r)
\nonumber\\
& =\max\left(\sqrt{p},\sqrt{p}\cdot\sqrt{\frac{\eta\left(\SU_{o}\right)r}{n}}\right)
\qquad \mbox{(see \eqref{eq:pre sigma Ais leg 1} and \eqref{eq:pre sigma Ais leg 2})}
\nonumber\\
&
=\sqrt{p}.\qquad\left(\eta\left(\SU_{o}\right)\le\frac{n}{r}\right)
\label{eq:sigma Ais 2}
\end{align}
We now compute
\begin{align}
\max\left(\log n\cdot \beta,\sqrt{\log n}\cdot\sigma\right)
& =\max\left(\log n\sqrt{\frac{\eta\left(\SU_{o}\right)r}{n}},\sqrt{{\log n\cdot p}}\right)
\qquad \mbox{(see \eqref{eq:b Ais 2} and  \eqref{eq:sigma Ais 2})}
\nonumber\\
& =\sqrt{{\log n\cdot p}}.\qquad\left(\mbox{if }p\ge \frac{\log n\cdot\eta\left(\SU_{o}\right)r}{n}\right)
\label{eq:other m cnd}
\end{align}
Therefore, for $\alpha\ge1$ and except with a probability of at most
$e^{-\alpha}$, it holds that
\begin{align}
\left\Vert S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}-\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}\right]\right\Vert &
\lesssim\alpha\max\left(\log n\cdot \beta,\sqrt{\log n}\cdot\sigma\right)
\qquad \mbox{(see Lemma \ref{lem:Bernie for spec})}
\nonumber\\
& =\alpha\sqrt{{\log n\cdot p}}.\label{eq:large dev 2}
\end{align}
Overall, substituting the large deviation bounds (\ref{eq:large dev 1})
and (\ref{eq:large dev 2}) into (\ref{eq:key}), we find that
\begin{align*}
p\left\Vert \widehat{Z}_{j}\right\Vert  & \le\left\Vert S_{o}^{*}P_{\omega_{j}}S_{o}-\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}\right]\right\Vert \cdot\left\Vert \widehat{Z}_{j}\right\Vert +\left\Vert S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}-\mathbb{E}\left[S_{o}^{*}P_{\omega_{j}}S_{o}^{\perp}\right]\right\Vert
\qquad \mbox{(see \eqref{eq:key})}
\\
 & \lesssim \alpha \sqrt{\frac{{\log r\cdot p \cdot \eta\left(\SU_{o}\right)r}}{n}}\cdot\left\Vert \widehat{Z}_{j}\right\Vert +\alpha\sqrt{{\log n\cdot p}},
 \qquad \mbox{(see \eqref{eq:large dev 1} and \eqref{eq:large dev 2})}
\end{align*}
except with a probability of at most $2e^{-\alpha}$ and under \eqref{eq:max b sigma Ais} and \eqref{eq:other m cnd}. It immediately follows that
\begin{align}
 \left\Vert \widehat{Z}_{j}\right\Vert
& \lesssim\frac{\alpha\sqrt{\frac{\log n}{ p}}}{1-\sqrt{\frac{\alpha^{2}\log r\cdot\eta\left(\SU_{o}\right)r}{pn}}}
\qquad \mbox{(see the next line)}
\nonumber\\
& \lesssim\alpha\sqrt{\frac{\log n}{p}},\qquad\left(\mbox{if }\frac{\alpha^{2}\log r\cdot\eta\left(\SU_{o}\right)r}{pn}\lesssim 1\right)
\end{align}
except with a probability of at most $2e^{-\alpha}$.
 In light of \eqref{eq:max b sigma Ais} and \eqref{eq:other m cnd}, we assume that $p \gtrsim \alpha^2 \log n \cdot \eta(\SU_o) r/n$.
Then, using
the union bound and with the choice of $\alpha=\alpha'\log b$, it
follows that
\[
\max_{j\in [1:b]}\left\Vert \widehat{Z}_{j}\right\Vert \lesssim\alpha'\log b\sqrt{\frac{\log n }{p}},
\]
provided that $p\gtrsim\alpha'^{2}\log^{2}b\cdot\log n\cdot\eta(\SU_{o})r/n$
and except with a probability of at most $2be^{-\alpha'\log b}=2 b^{1-\alpha'}$. Finally,
invoking (\ref{eq:brk .51}), we conclude that
\begin{align*}
\left\Vert O\left(P_{\SU_{o}^{\perp}}S_{Q}\right)\right\Vert _{F} &
\le \max_j \l\| \widehat{Z}_j \r\| \cdot
\l\| P_{\SU_o^\perp} P_{\SU} \r\|_F \|Q\|
\qquad \mbox{(see \eqref{eq:brk .51})}
\nonumber\\
& \lesssim\alpha'\log b\sqrt{\frac{\log n}{p}}\cdot\left\Vert P_{\SU_{o}^{\perp}}P_{\SU}\right\Vert _{F}\|Q\|,
\end{align*}
which completes the proof of Lemma \ref{lem:O norm}.


\section{Properties of a Standard Random Gaussian Matrix}
\label{sec:cnd coh of Gaussian}

As a supplement to Remark \ref{rem:coeff}, we show here that a standard  random Gaussian matrix $G\in\mathbb{R}^{b\times r}$ is well-conditioned and incoherent, when $b \gtrsim r$. Let $\sigma_{\max}(G)\ge \sigma_{\min}(G)$ denote the largest and smallest singular values of $G$, respectively.  From \cite[Corollary 5.35]{vershynin2010introduction} and for fixed $\alpha \ge 1$, recall that
\begin{align}
\label{eq:norms of Gaussian}
\sqrt{b}- \Cl{Gaussian} \alpha \sqrt{r} \le \sigma_{\min}(G)\le \sigma_{\max}(G)\le \sqrt{b} + \Cr{Gaussian} \alpha \sqrt{r},
\end{align}
for an absolute constant $\Cr{Gaussian}>0$ and except with a probability of at most $e^{-\alpha^2 r}$. It follows that
\begin{align}
\nu(G) = \frac{\sigma_{\max}(G)}{\sigma_{\min}(G)}
\le \frac{\sqrt{b}+\Cr{Gaussian}\alpha \sqrt{r}}{\sqrt{b}-\Cr{Gaussian}\alpha \sqrt{r}},
\end{align}
which, if $\sqrt{b}\ge 2\Cr{Gaussian}\alpha \sqrt{r}$, reduces to $\nu(G)=O(1)$ with high probability.

For the coherence, note that $G(G^*G)^{-\frac{1}{2}}\in\mathbb{R}^{b\times r}$ is an orthonormal basis for $\mbox{span}(G)$. Using the definition of coherence, we write that
\begin{align}
\eta\l( \mbox{span}(G) \r) & = \frac{b}{r} \max_{i\in[1:b] }
\l\| G[i,:] \l(G^*G \r)^{-\frac{1}{2}}\r\|_2^2
\qquad \mbox{(see \eqref{eq:def of coh})} \nonumber\\
& \le \frac{b}{r} \max_i \l\| G[i,:]\r\|_2^2 \cdot
\l\|\l(G^*G\r)^{-\frac{1}{2}} \r\|^2 \nonumber\\
& = \frac{b}{r} \max_i \l\| G[i,:]\r\|_2^2 \cdot \l( \sigma_{\min}(G)\r)^{-2} \nonumber\\
& \le \frac{b}{r}\max_i \l\| G[i,:]\r\|_2^2 \cdot \l(\sqrt{b}-\Cr{Gaussian}\alpha \sqrt{r}\r)^{-2}
\qquad \mbox{(see \eqref{eq:norms of Gaussian})} \nonumber\\
& \le \frac{b}{r} \max_i \l\| G[i,:]\r\|_2^2 \cdot \l(\frac{b}{2}-\Cr{Gaussian}^2\alpha^2 r \r)^{-1}
\qquad \l( (a-b)^2 \ge \frac{a^2}{2}-b^2\r) \nonumber\\
& \lesssim  \frac{b}{r} \max_i \l\| G[i,:]\r\|_2^2 \cdot \l(b-\Cl{Gaussian2} \alpha^2 r \r)^{-1},
\label{eq:eta Gaussian 1}
\end{align}
for a constant $\Cr{Gaussian2}>0$ and
except with a probability of at most $e^{-\alpha^2 r}$. For fixed $i$, $\|G[i,:]\|_2^2$ is a chi-squared random variable  with $r$ degrees of freedom so that
\begin{align}
\Pr\l[\l\| G[i,:]\r\|_2^2 \gtrsim t \cdot r \r] \le e^{-t},
\end{align}
for any $t\ge 1$. An application of the union bound and the choice of  $t=\Cl{GaussianCoh} \alpha \log b$ then lead us to
\begin{align}
\Pr\l[\max_{i\in[1:b]}\l\| G[i,:]\r\|_2^2 \gtrsim \Cr{GaussianCoh} \alpha \log b\cdot r \r] \le b^{-(\Cr{GaussianCoh}\alpha -1)}.
\label{eq:union bnd Gaussian}
\end{align}
Substituting the bound above back into \eqref{eq:eta Gaussian 1} yields that
\begin{align}
\eta\l( \mbox{span}(G) \r)  & \lesssim \frac{b}{r}\max_i \l\| G[i,:]\r\|_2^2 \cdot \l(b-\Cr{Gaussian2}r\r)^{-1} \qquad \mbox{(see \eqref{eq:eta Gaussian 1})} \nonumber\\
& \lesssim  \frac{b}{r} \cdot  r\alpha \log b\cdot \l(b-\Cr{Gaussian2}r\r)^{-1}
\qquad \mbox{(see \eqref{eq:union bnd Gaussian})} \nonumber\\
& = \frac{\alpha b\log b}{b-\Cr{Gaussian2}\alpha^2 r},
\end{align}
except with a probability of at most $e^{-\alpha r}+b^{-(\Cr{GaussianCoh}\alpha -1)}$. In particular, when $b \ge  2  \Cr{Gaussian2} \alpha^2 r$, we find that $\eta(\mbox{span}(G))\lesssim \alpha \log b$ with high probability.

\end{document}

SNIPE_v10.bbl
\begin{thebibliography}{10}

\bibitem{van2012subspace}
P.~van Overschee and B.~L. de~Moor.
\newblock {\em Subspace identification for linear systems: {T}heory,
  implementation, applications}.
\newblock Springer US, 2012.

\bibitem{ardekani1999activation}
B.~A. Ardekani, J.~Kershaw, K.~Kashikura, and I.~Kanno.
\newblock Activation detection in functional {MRI} using subspace modeling and
  maximum likelihood estimation.
\newblock {\em IEEE Transactions on Medical Imaging}, 18(2):101--114, 1999.

\bibitem{krim1996two}
H.~Krim and M.~Viberg.
\newblock Two decades of array signal processing research: {T}he parametric
  approach.
\newblock {\em IEEE Signal processing magazine}, 13(4):67--94, 1996.

\bibitem{tong1998multichannel}
L.~Tong and S.~Perreau.
\newblock Multichannel blind identification: {F}rom subspace to maximum
  likelihood methods.
\newblock {\em Proceedings of IEEE}, 86:1951--1968, 1998.

\bibitem{cunningham2015linear}
J.~P. Cunningham and Z.~Ghahramani.
\newblock Linear dimensionality reduction: Survey, insights, and
  generalizations.
\newblock {\em Journal of Machine Learning Research}, 16:2859--2900, 2015.

\bibitem{recht2011hogwild}
B.~Recht, C.~Re, S.~J. Wright, and F.~Niu.
\newblock Hogwild: {A} lock-free approach to parallelizing stochastic gradient
  descent.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
  693--701, 2011.

\bibitem{nocedal2006numerical}
J.~Nocedal and S.~J. Wright.
\newblock {\em Numerical optimization}.
\newblock Springer Series in Operations Research and Financial Engineering.
  Springer New York, 2006.

\bibitem{golub2013matrix}
G.~H. Golub and C.~F. Van~Loan.
\newblock {\em Matrix computations}.
\newblock Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins
  University Press, 2013.

\bibitem{balzano2015local}
L.~Balzano and S.~J. Wright.
\newblock Local convergence of an algorithm for subspace identification from
  partial data.
\newblock {\em Foundations of Computational Mathematics}, 15(5):1279--1314,
  2015.

\bibitem{mitliagkas2014streaming}
I.~Mitliagkas, C.~Caramanis, and P.~Jain.
\newblock Streaming {PCA} with many missing entries.
\newblock {\em Preprint}, 2014.

\bibitem{chen2015incoherence}
Y.~Chen.
\newblock Incoherence-optimal matrix completion.
\newblock {\em IEEE Transactions on Information Theory}, 61(5):2909--2923,
  2015.

\bibitem{olver2010nist}
F.~W.~J. Olver, National~Institute of~Standards, and Technology (U.S.).
\newblock {\em NIST Handbook of Mathematical Functions Hardback and CD-ROM}.
\newblock Cambridge University Press, 2010.

\bibitem{balzano2010online}
L.~Balzano, R.~Nowak, and B.~Recht.
\newblock Online identification and tracking of subspaces from highly
  incomplete information.
\newblock In {\em Annual Allerton Conference on Communication, Control, and
  Computing (Allerton)}, pages 704--711. IEEE, 2010.

\bibitem{halko2011finding}
N.~Halko and J.~A. Martinsson, P. G.and~Tropp.
\newblock Finding structure with randomness: {P}robabilistic algorithms for
  constructing approximate matrix decompositions.
\newblock {\em SIAM review}, 53(2):217--288, 2011.

\bibitem{ge2016matrix}
R.~Ge, J.~D. Lee, and T.~Ma.
\newblock Matrix completion has no spurious local minimum.
\newblock {\em arXiv preprint arXiv:1605.07272}, 2016.

\bibitem{chi2013petrels}
Y.~Chi, Y.~C. Eldar, and R.~Calderbank.
\newblock {PETRELS}: {P}arallel subspace estimation and tracking by recursive
  least squares from partial observations.
\newblock {\em IEEE Transactions on Signal Processing}, 61(23):5947--5959,
  2013.

\bibitem{mardani2015subspace}
M.~Mardani, G.~Mateos, and G.~B. Giannakis.
\newblock Subspace learning and imputation for streaming big data matrices and
  tensors.
\newblock {\em IEEE Transactions on Signal Processing}, 63(10):2663--2677,
  2015.

\bibitem{xie2013change}
Y.~Xie, J.~Huang, and R.~Willett.
\newblock Change-point detection for high-dimensional time series with missing
  data.
\newblock {\em IEEE Journal of Selected Topics in Signal Processing},
  7(1):12--27, 2013.

\bibitem{balzano2013grouse}
L.~Balzano and S.~J Wright.
\newblock On {GROUSE} and incremental {SVD}.
\newblock In {\em IEEE International Workshop on Computational Advances in
  Multi-Sensor Adaptive Processing (CAMSAP)}, pages 1--4. IEEE, 2013.

\bibitem{bunch1978updating}
J.~R. Bunch and C.~P. Nielsen.
\newblock Updating the singular value decomposition.
\newblock {\em Numerische Mathematik}, 31(2):111--129, 1978.

\bibitem{balsubramani2013fast}
A.~Balsubramani, S.~Dasgupta, and Y.~Freund.
\newblock The fast convergence of incremental {pca}.
\newblock In {\em Advances in Neural Information Processing Systems}, pages
  3174--3182, 2013.

\bibitem{oja1985stochastic}
E.~Oja and J.~Karhunen.
\newblock On stochastic approximation of the eigenvectors and eigenvalues of
  the expectation of a random matrix.
\newblock {\em Journal of mathematical analysis and applications},
  106(1):69--84, 1985.

\bibitem{lounici2014high}
K.~Lounici.
\newblock High-dimensional covariance matrix estimation with missing
  observations.
\newblock {\em Bernoulli}, 20(3):1029--1058, 2014.

\bibitem{candes2009exact}
E.~J. Cand{\`e}s and B.~Recht.
\newblock Exact matrix completion via convex optimization.
\newblock {\em Foundations of Computational mathematics}, 9(6):717--772, 2009.

\bibitem{eftekhari2016mc}
A.~Eftekhari, M.~B. Wakin, and R.~A. Ward.
\newblock Mc2: {A} two-phase algorithm for leveraged matrix completion.
\newblock {\em arXiv preprint arXiv:1609.01795}, 2016.

\bibitem{gross2011recovering}
D.~Gross.
\newblock Recovering low-rank matrices from few coefficients in any basis.
\newblock {\em IEEE Transactions on Information Theory}, 57(3):1548--1566,
  2011.

\bibitem{tropp2012user}
J.~A. Tropp.
\newblock User-friendly tail bounds for sums of random matrices.
\newblock {\em Foundations of Computational Mathematics}, 12(4):389--434, 2012.

\bibitem{wedin1972perturbation}
P.~Wedin.
\newblock Perturbation bounds in connection with singular value decomposition.
\newblock {\em BIT Numerical Mathematics}, 12(1):99--111, 1972.

\bibitem{vershynin2010introduction}
R.~Vershynin.
\newblock Introduction to the non-asymptotic analysis of random matrices.
\newblock In Y.~C. Eldar and G.~Kutyniok, editors, {\em Compressed Sensing:
  Theory and Applications}, pages 95--110. Cambridge University Press, 2012.

\end{thebibliography}
