text:
\section{Introduction}

Least squares approximations of the form
\[
\min_{{\mathbf{x}}\in\mathbb{R}^{d}}\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}-{\mathbf{b}}}
\]
are fundamental building blocks in computational science and statistical
data analysis, with applications ranging from statistical data analysis
to inverse problems. However, it is well appreciated, especially in
the aforementioned application areas, that regularization is often
the key to achieving the best results.

One of the basic methods for regularizing least squares approximations
is principal component regression (PCR)~\cite{hotelling1933analysis,kendall1957course,artemiou2009principal}.
Given a data matrix $\ensuremath{{\bm{\mathrm{A}}}}$, a right hand side ${\mathbf{b}}$ and a target
rank $k$, PCR is computed by first computing the coefficients $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$
corresponding to the top $k$ principal components of $\ensuremath{{\bm{\mathrm{A}}}}$ (i.e.,
to dominant right invariant subspace of $\ensuremath{{\bm{\mathrm{A}}}}$), then regressing
on $\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$ and ${\mathbf{b}}$, and finally projecting the solution
back to the original space. In short, the PCR estimator is ${\mathbf{x}}_{k}=\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})^{\textsc{+}}{\mathbf{b}}$
and regularization is achieved via PCA based dimensionality reduction.
While there is some criticism of PCR in the statistical literature~\cite{artemiou2009principal,jolliffe1982note},
it is nevertheless a valuable tool in the toolbox of practitioners.

Up until recent breakthroughs on fast methods for least squares approximations,
there was little penalty in terms of computational complexity when
switching from ordinary least squares (OLS) to PCR. Indeed, the complexity
of SVD based computation of the dominant invariant subspace is $O(nd\min(n,d))$,
and this matches the asymptotic complexity of straightforward computation
of the OLS solution (i.e., via direct methods). However, recent progress
on fast sketching based algorithms for linear regression~\cite{DMMS11,RT08,MSM14,CW17,woodruff2014sketching}
has created a gap: exact computation of the principal components still
requires SVD so the overall complexity is still $O(nd\min(n,d))$,
even though the OLS stage is faster. The gap is not insubstantial:
when learning with large scale data (either large $n$, or large $d$),
$O(nd\min(n,d))$ is often infeasible, but modern sketching based
linear regression methods are.

\subsection{Contributions}

In this paper, we study the use of dimensionality reduction prior
to computing PCR (so we can compute PCR on a smaller input matrix).
In particular, for a data matrix $\ensuremath{{\bm{\mathrm{A}}}}$, we relate the PCR solution
of $\ensuremath{{\bm{\mathrm{A}}}}\mat R$, where $\mat R$ is any dimensionality reduction matrix,
to the PCR solution of $\ensuremath{{\bm{\mathrm{A}}}}$. To do so, we study the notion of
approximate PCR both from an optimization perspective and from a statistical
perspective, and provide conditions on $\mat R$ that guarantee that
after projecting the solution back to the full space (by multiplying
by $\mat R^{\textsc{T}}$) we have an approximate PCR solution with rigorous
statistical risk bounds. These results are described in Section~\ref{sec:dim-reduce-pcr}.

We then leverage the aforementioned results to design fast, sketching
based, algorithms for approximate PCR. We propose algorithms specialized
for the several cases (in the following, $n$ is number of data points,
$d$ is dimension of the data): large $n$ (using left sketching),
large $d$ (using right sketching), and both $n$ and $d$ large (using
two-sided sketching). Furthermore, we propose an input-sparsity time
algorithm for approximate PCR. These results are described in Section~\ref{sec:alg}.

We also consider computing approximate PCR in the streaming model,
providing the first algorithm for computing approximate PCR in a stream.
We also provide a fast algorithm for approximate Kernel PCR (polynomial
kernel only). These results are described in Section~\ref{sec:extensions}.

Finally, empirical results (Section \ref{sec:experiments}) clearly
demonstrate the ability of our proposed algorithms to compute approximate
PCR solution, the correctness of our theoretical analysis, and the
advantages of using our techniques instead of simpler techniques like
compressed least squares.

In general, unlike previous works on randomized methods for PCR (which
we discuss in the next subsection), we analyze the use of sketching
for PCR from a sketch-and-solve approach. We discuss the various advantages
and disadvantages of the sketch-and-solve approach in comparison to
iterative based approaches, in the next subsection.

\subsection{Related Work}

Recently matrix sketching, such as the use of random projections,
has emerged as a powerful technique for accelerating and scaling many
important statistical learning techniques. See recent surveys by Woodruff~\cite{woodruff2014sketching}
and Mahoney et al.~\cite{YMM15} for an extensive exposition on this
subject. So far, there has been limited research on the use of matrix
sketching in the context of principal component regression.

One natural strategy for leveraging sketching in the context of PCR
is to use approximate principal components. Approximate principal
components can be computed using fast sketching based algorithm for
approximate PCA (also known as 'randomized SVD')~\cite{halko2011finding,woodruff2014sketching}.
This was recently explored by Boutsidis and Magdon-Ismail~\cite{BM14}.
The authors show that the if the number of subspace iterations is
sufficiently large, one obtain a bound on the sub-optimality of the
approximate solution and on the error of the solution vector. We too
bound the sub-optimality of our solutions, but instead of bounding
the error of the solution vector, we bound their distance to the right
dominant subspace, or bound the distance of the projection to the
left dominant subspace.

Frostig et al. leverage fast randomized algorithms for ridge regression
to design iterative algorithms for principal component regression
and principal component projection~\cite{frostig2016principal}.
Forstig et al.'s results were later improved by Allen-Zhu and Li~\cite{AL17}.
Both of the aforementioned methods use iterations, while our work
explores the use of a sketch-and-solve approach. While it is true
that better accuracies can be achieved using iterative methods with
sketching based accelerators~\cite{RT08,AMT10,MSM14,GOS16,ACW17},
there are some advantages in using a sketch-and-solve approach. In
particular, sketch-and-solve algorithms are typically faster. However
this comes at the cost: sketch-and-solve algorithms typically provide
cruder approximations. Nevertheless, it is not uncommon for these
cruder approximations to be sufficient in machine learning applications.
Another advantage of the sketch-and-solve approach is that it is more
amenable to streaming and kernelization; we consider both in this
paper.

Closely related to our work is recent work on Compressed Least Squares
(CLS)~\cite{maillard2009compressed,Kaban14,slawski2017compressed,Slawski17,THM17}.
In particular, our statistical analysis (section~\ref{subsec:statistical})
is inspired by recent statistical analysis of CLS~\cite{slawski2017compressed,Slawski17,Kaban14}.
Additionally, CLS is sometimes considered as a computationally attractive
alternative to PCR~\cite{Slawski17,THM17}. While CLS certainly uses
matrix sketching to compress the matrix, it also uses the compression
to regularize the problem. The mix between compression for scalability
and compression for regularization reduces the ability to fine tune
the method to the needs at hand, and thus obtain the best possible
results. In contrast, our methods uses sketching primarily to approximate
the principal components and as such serves as a means for scalability
only. We propose methods that are computationally as attractive as
CLS, and are more faithful to the behavior of PCR (in fact, CLS is
a special case of one of our proposed algorithms). These advantages
over CLS are also evident in the experimental results reported in
Section~\ref{sec:experiments}.

Principal component regression is a form of least squares regression
with convex constraints (once the dominant subspace has been found).
Pilanchi and Wainwright recently explored the effect of regularization
on the sketch size for least squares regression~\cite{pilanci2015randomized,pilanci2016iterative}.
In the aforementioned papers, sketching is applied only to the objective,
while the constraint is enforced exactly. This is unsatisfactory in
the context of PCR since for PCR the constraints are gleaned from
the input, and enforcing them is as expensive as solving the problem
exactly. In contrast, our method uses sketching not only to compress
the objective function, but also to approximate the constraint set.

Ridge regression (also known as Tikhonov regularization) is another
popular and well studied method for regularizing least squares solutions.
It also closely related to PCR in the sense that the ridge term can
be viewed as a soft damping of the singular values. Recently several
sketching-based algorithms have been suggested to accelerate the solution
of ridge regression~\cite{ChenEtAl15,ACW17b,WGM17,CYD18}.

\section{Preliminaries}

\subsection{Notation and Basic Definitions}

We denote scalars using Greek letters or using $x,y,\dots$. Vectors
are denoted by ${\mathbf{x}},{\mathbf{y}},\dots$ and matrices by $\ensuremath{{\bm{\mathrm{A}}}},\mat B,\dots$.
The $s\times s$ identity matrix is denoted $\mat I_{s}$. We use the
convention that vectors are column-vectors. $\nnz{\ensuremath{{\bm{\mathrm{A}}}}}$ denotes
the number of non-zeros in $\ensuremath{{\bm{\mathrm{A}}}}$. The notation $\alpha=(1\pm\gamma)\beta$
means that $(1-\gamma)\beta\leq\alpha\leq(1+\gamma)\beta$, and the
notation $\alpha=\beta\pm\gamma$ means that $|\alpha-\beta|\leq\gamma$.

Given a matrix $\ensuremath{{\bm{\mathrm{X}}}}\in\mathbb{R}^{m\times n}$, let $\ensuremath{{\bm{\mathrm{X}}}}=\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}\Sigma_{\ensuremath{{\bm{\mathrm{X}}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}^{\textsc{T}}$
be a \emph{thin SVD} of $\ensuremath{{\bm{\mathrm{X}}}}$, i.e. $\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}\in\mathbb{R}^{m\times\min(m,n)}$
is a matrix with orthonormal columns, $\Sigma_{\ensuremath{{\bm{\mathrm{X}}}}}\in\mathbb{R}^{\min(m,n)\times\min(m,n)}$
is a diagonal matrix with the non-negative singular values on the
diagonal, and $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}\in\mathbb{R}^{n\times\min(m,n)}$ is a matrix
with orthonormal columns. The thin SVD decomposition is not necessarily
unique, so when we use this notation we mean that the statement is
correct for any such decomposition. A thin SVD decomposition can be
computed in $O(mn\min(m,n))$. We denote the singular values of $\ensuremath{{\bm{\mathrm{X}}}}$
by $\sigma_{\max}(\ensuremath{{\bm{\mathrm{X}}}})=\sigma_{1}(\ensuremath{{\bm{\mathrm{X}}}})\geq\dots\geq\sigma_{\min(m,n)}(\ensuremath{{\bm{\mathrm{X}}}})=\sigma_{\min}(\ensuremath{{\bm{\mathrm{X}}}})$,
omitting the matrix from the notation if the relevant matrix is clear
from the context. For $k\leq\min(m,n)$, we use $\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}},k}$
(respectively $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}},k})$ to denote the matrix consisting
of the first $k$ columns of $\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}$ (respectively $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}$),
and use $\Sigma_{\ensuremath{{\bm{\mathrm{X}}}},k}$ to denote the leading $k\times k$ minor
of $\Sigma_{\ensuremath{{\bm{\mathrm{X}}}}}$. We use $\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}},k+}$ (respectively $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}},k+})$
to denote the matrix consisting of the last $\min(m,n)-k$ columns
of $\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}$ (respectively $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}$), and use $\Sigma_{\ensuremath{{\bm{\mathrm{X}}}},k+}$
to denote the lower-right $(\min(m,n)-k)\times(\min(m,n)-k)$ block
of $\Sigma_{\ensuremath{{\bm{\mathrm{X}}}}}$. In other words,
\[
\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}=\left[\begin{array}{cc}
\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}},k} & \ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}},k+}\end{array}\right]\quad\Sigma_{\ensuremath{{\bm{\mathrm{X}}}}}=\left[\begin{array}{cc}
\Sigma_{\ensuremath{{\bm{\mathrm{X}}}},k} & 0\\
0 & \Sigma_{\ensuremath{{\bm{\mathrm{X}}}},k+}
\end{array}\right]\quad\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}=\left[\begin{array}{cc}
\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}},k} & \ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}},k+}\end{array}\right].
\]
The \emph{Moore-Penrose pseudo-inverse} of $\ensuremath{{\bm{\mathrm{X}}}}$ is $\ensuremath{{\bm{\mathrm{X}}}}^{\textsc{+}}\coloneqq\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}\Sigma_{\ensuremath{{\bm{\mathrm{X}}}}}^{\textsc{+}}\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}^{\textsc{T}}$
where $\Sigma_{\ensuremath{{\bm{\mathrm{X}}}}}^{\textsc{+}}=\diag{\sigma_{1}(\ensuremath{{\bm{\mathrm{X}}}})^{\textsc{+}},\dots,\sigma_{\min(m,n)}(\ensuremath{{\bm{\mathrm{X}}}})^{\textsc{+}}}$
with $a^{+}=a^{-1}$ when $a\neq0$ and $0$ otherwise.

The \emph{stable rank }of a matrix $\ensuremath{{\bm{\mathrm{X}}}}$ is $\sr{\ensuremath{{\bm{\mathrm{X}}}}}\coloneqq\FNormS{\ensuremath{{\bm{\mathrm{X}}}}}/\TNormS{\ensuremath{{\bm{\mathrm{X}}}}}$.
The $k$-th \emph{relative gap} of a matrix $\ensuremath{{\bm{\mathrm{X}}}}$ is
\[
\gap{\ensuremath{{\bm{\mathrm{X}}}}}k=\frac{\sigma_{k}^{2}-\sigma_{k+1}^{2}}{\sigma_{1}^{2}}\,.
\]

For a subspace ${\cal U}$, we use $\ensuremath{{\bm{\mathrm{P}}}}_{{\cal U}}$ to denote the
\emph{orthogonal projection matrix} onto ${\cal U}$, and $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}$
for the projection matrix on the column space of $\ensuremath{{\bm{\mathrm{X}}}}$ (i.e. $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}=\ensuremath{{\bm{\mathrm{P}}}}_{\range{\ensuremath{{\bm{\mathrm{X}}}}}}$).
We have $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}=\ensuremath{{\bm{\mathrm{X}}}}\matX^{\textsc{+}}$. The \emph{complementary
projection matrix }is $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}^{\perp}=\mat I-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{X}}}}}$.
A useful property of projection matrices is that if ${\cal S}\subseteq{\cal T}$
then $\ensuremath{{\bm{\mathrm{P}}}}_{{\cal S}}\ensuremath{{\bm{\mathrm{P}}}}_{{\cal T}}=\ensuremath{{\bm{\mathrm{P}}}}_{{\cal T}}\ensuremath{{\bm{\mathrm{P}}}}_{{\cal S}}=\ensuremath{{\bm{\mathrm{P}}}}_{{\cal S}}$.
Furthermore, we note the following result.
\begin{thm}
[Theorem 2.3 in \cite{Stewart_perturbation}]\label{thm:projection_equality}
For any $\ensuremath{{\bm{\mathrm{A}}}}$ and $\ensuremath{{\bm{\mathrm{B}}}}$ with the same number of rows, the following
statements hold:
\begin{enumerate}
\item If $\rank{\ensuremath{{\bm{\mathrm{A}}}}}=\rank{\ensuremath{{\bm{\mathrm{B}}}}}$, then the singular values of $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{B}}}}}^{\perp}$
and $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{B}}}}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}^{\perp}$ are the same, so
\[
\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{B}}}}}^{\perp}}=\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{B}}}}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}^{\perp}}
\]
\item Moreover the nonzero singular values $\sigma$ of $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{B}}}}}^{\perp}$
correspond to pairs $\pm\sigma$ of eigenvalues of $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{B}}}}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}$,
so
\[
\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{B}}}}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}}=\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{B}}}}}^{\perp}}
\]
\item If $\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{B}}}}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}}<1$, then $\rank{\ensuremath{{\bm{\mathrm{A}}}}}=\rank{\ensuremath{{\bm{\mathrm{B}}}}}$.
\end{enumerate}
\end{thm}

\subsection{Principal Component Regression and Principal Component Projection }

In the \emph{Principal Component Regression (PCR)} problem, we are
given an input $n\textrm{-by-}d$ data matrix $\ensuremath{{\bm{\mathrm{A}}}}$, a right hand
side ${\mathbf{b}}\in\mathbb{R}^{n}$, and a rank parameter $k$ which is smaller or
equal to the rank of $\ensuremath{{\bm{\mathrm{A}}}}$. Furthermore, we assume that there is
an non-zero eigengap at $k$: $\sigma_{k}>\sigma_{k+1}$. The goal
is to find the PCR solution, ${\mathbf{x}}_{k}$, defined as
\begin{equation}
{\mathbf{x}}_{k}\coloneqq\arg\min_{{\mathbf{x}}\in\range{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}}\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}-{\mathbf{b}}}.\label{eq:pcr-opt-problem}
\end{equation}
It is easy to verify that ${\mathbf{x}}_{k}=\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})^{\textsc{+}}{\mathbf{b}}=\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}\Sigma_{\ensuremath{{\bm{\mathrm{A}}}},k}^{-1}\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\textsc{T}}{\mathbf{b}}$.
The \emph{Principal Component Projection (PCP) }of ${\mathbf{b}}$ is ${\mathbf{b}}_{k}\coloneqq\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}=\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}{\mathbf{b}}$.

Straightforward computation of ${\mathbf{x}}_{k}$ and ${\mathbf{b}}_{k}$ via the SVD
takes $O(nd\min(n,d))$ operations\footnote{The complexity when using iterative algorithms (e.g. Lanczos) to compute
only the dominant invariant spaces depend on several additional facts
and in particular on spectral properties of the matrix and sparsity
level. Thus, to avoid overly complicating the discussion on computational
complexity, we refrain from further discussion of iterative methods
for computing dominant eigenspaces}. We are primarily interested in finding faster algorithms that compute
an approximate PCR or PCP solution (we formalize the terms 'approximate
PCP/PCR' in Section~\ref{sec:dim-reduce-pcr}). Throughout the paper,
we use $\ensuremath{{\bm{\mathrm{A}}}},{\mathbf{b}},$ and $k$ as the arguments of the PCR/PCP problem
to be solved.

\subsection{\label{subsec:Matrix-Perturbations}Matrix Perturbations and Distance
Between Subspaces}

Our analysis uses matrix perturbation theory extensively. We now describe
the basics of this theory and the results we use.

The \emph{principal angles} $\theta_{j}\in[0,\pi/2]$ between two
subspaces ${\cal U}$ and ${\cal W}$ are recursively defined by the
identity
\[
\cos(\theta_{j})=\max_{{\mathbf{u}}\in{\cal U}}\max_{{\mathbf{w}}\in{\cal W}}{\mathbf{u}}^{\textsc{T}}{\mathbf{w}}\,\text{s.t.}\,\TNorm{{\mathbf{u}}}=1,\TNorm{{\mathbf{w}}}=1,\forall i<j.{\mathbf{u}}_{i}^{\textsc{T}}{\mathbf{u}}=0,{\mathbf{w}}_{i}^{\textsc{T}}{\mathbf{w}}=0\,.
\]
We use ${\mathbf{u}}_{j}$ and ${\mathbf{w}}_{j}$ to denote the vectors for which $\cos(\theta_{j})={\mathbf{u}}_{j}^{\textsc{T}}{\mathbf{w}}_{j}$.
Let $\Theta({\cal {\cal U}},{\cal W})$ denote the $d\times d$ diagonal
matrix whose $j$th diagonal entry is the $j$th principal angle,
and as usual we allow writing matrices instead of subspaces as short-hand
for the column space of the matrix. Henceforth, when we write a function
on $\Theta(\cdot,\cdot)$, i.e. $\sin(\Theta(\ensuremath{{\bm{\mathrm{U}}}},\text{\ensuremath{\mat W}}$)),
we mean evaluating the function entrywise on the diagonal only. It
is well known~ \cite[section 6.4.3]{golub2012matrix} that if $\ensuremath{{\bm{\mathrm{U}}}}$
(respectively $\mat W$) is a matrix with orthonormal columns whose
column space is equal to ${\cal U}$ (respectively ${\cal W}$) then
\[
\sigma_{j}(\ensuremath{{\bm{\mathrm{U}}}}^{\textsc{T}}\mat W)=\cos(\theta_{j}).
\]
The following lemma connects the tangent of the principal angles to
the spectral norm of an appropriate matrix.
\begin{lem}
[Lemma 4.3 in \cite{drineas2016structural}] \label{lem:tan-to-spectral}Let
$\ensuremath{{\bm{\mathrm{Q}}}}\in\mathbb{R}^{n\times s}$ have orthonormal columns, and let $\mat W=(\begin{array}{cc}
\mat W_{k} & \mat W_{k+}\end{array})\in\mathbb{R}^{n\times n}$ be an orthogonal matrix where $\mat W_{k}\in\mathbb{R}^{n\times k}$ with
$k\leq s$. If $\rank{\mat W_{k}^{\textsc{T}}\ensuremath{{\bm{\mathrm{Q}}}}}=k$ then
\[
\TNorm{\tan\Theta(\ensuremath{{\bm{\mathrm{Q}}}},\mat W_{k})}=\TNorm{(\mat W_{k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{Q}}}})(\mat W_{k}^{\textsc{T}}\ensuremath{{\bm{\mathrm{Q}}}})}\,.
\]
\end{lem}

Matrix perturbation theory studies how a perturbation of a matrix
translate to perturbations of the matrix's eignevalues and eigenspaces.
In order to bound the perturbation of an eigenspace, one needs some
notion of distance between two subspaces. One common distance metric
between two subspaces is
\begin{equation}
d_{2}({\cal U},{\cal W})\coloneqq\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{{\cal U}}-\ensuremath{{\bm{\mathrm{P}}}}_{{\cal W}}}\,.\label{eq:distance_between_subspaces}
\end{equation}
If $\ensuremath{{\bm{\mathrm{U}}}}$ and $\ensuremath{{\bm{\mathrm{V}}}}$ have the same number of columns, and both
have orthonormal columns, then
\[
d_{2}(\ensuremath{{\bm{\mathrm{U}}}},\ensuremath{{\bm{\mathrm{V}}}})=\sqrt{1-\sigma_{\min}(\ensuremath{{\bm{\mathrm{U}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}})^{2}}=\sin(\theta_{\max})=\TNorm{\sin\Theta(\ensuremath{{\bm{\mathrm{U}}}},\ensuremath{{\bm{\mathrm{V}}}})}
\]
where $\theta_{\max}$ is the maximum principal angle between $\range{\ensuremath{{\bm{\mathrm{U}}}}}$
and $\range{\ensuremath{{\bm{\mathrm{V}}}}}$~\cite[section 6.4.3]{golub2012matrix}.

A classical result that bounds the distance between the dominant subspaces
of two symmetric matrices in terms of the spectral norm of difference
between the two matrices is the Davis-Kahan $\sin(\Theta)$ theorem~\cite[Section 2]{davis1970rotation}.
We need the following corollary of this theorem:
\begin{thm}
[Corollary of Davis-Kahan $\sin\Theta$ Theorem \cite{davis1970rotation}]\label{thm:sin-theta-1}
Let $\ensuremath{{\bm{\mathrm{A}}}},\tilde{\ensuremath{{\bm{\mathrm{A}}}}}\in\mathbb{R}^{n\times n}$ be two symmetric matrices,
both of rank at least $k$. Suppose that $\lambda_{k}>\tilde{\lambda}_{k+1}$
where $\lambda_{1}\geq\dots\geq\lambda_{n}$ and $\tilde{\lambda}_{1}\geq\dots\geq\tilde{\lambda}_{n}$
are the eigenvalues of $\ensuremath{{\bm{\mathrm{A}}}}$ and $\tilde{\ensuremath{{\bm{\mathrm{A}}}}}$. We have
\[
d_{2}(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k},\ensuremath{{\bm{\mathrm{V}}}}_{\tilde{\ensuremath{{\bm{\mathrm{A}}}}},k})\leq\frac{\TNorm{\ensuremath{{\bm{\mathrm{A}}}}-\tilde{\ensuremath{{\bm{\mathrm{A}}}}}}}{\lambda_{k}-\tilde{\lambda}_{k+1}}\,.
\]
\end{thm}

\begin{proof}
We use the following variant of the $\sin\Theta$ Theorem (see~\cite[Theorem 2.16]{stewart2001matrix}):
suppose a symmetric matrix $\ensuremath{{\bm{\mathrm{B}}}}$ has a spectral representation
\[
\ensuremath{{\bm{\mathrm{B}}}}=\ensuremath{{\bm{\mathrm{X}}}}\mat L\ensuremath{{\bm{\mathrm{X}}}}^{\textsc{T}}+\mat Y\ensuremath{{\bm{\mathrm{M}}}}\mat Y^{\textsc{T}}
\]
where $\left[\ensuremath{{\bm{\mathrm{X}}}}\,\mat Y\right]$ is square orthonormal. Let the
orthonormal matrix $\mat Z$ be of the same dimensions as $\ensuremath{{\bm{\mathrm{X}}}}$
and suppose that
\[
\mat R=\ensuremath{{\bm{\mathrm{B}}}}\mat Z-\mat Z\ensuremath{{\bm{\mathrm{N}}}}
\]
where $\ensuremath{{\bm{\mathrm{N}}}}$ is symmetric. Furthermore, suppose that the spectrum
of $\ensuremath{{\bm{\mathrm{N}}}}$ is contained in some interval $[\alpha,\beta]$ and that
for some $\delta>0$ the spectrum of $\ensuremath{{\bm{\mathrm{M}}}}$ lies outside of $[\alpha-\delta,\beta+\delta]$.
Then,
\[
\TNorm{\sin\Theta(\ensuremath{{\bm{\mathrm{X}}}},\mat Z)}\leq\frac{\TNorm{\mat R}}{\delta}\,.
\]
We prove Theorem~\ref{thm:sin-theta-1} by applying the aforementioned
variant of the $\sin\Theta$ Theorem with: $\ensuremath{{\bm{\mathrm{B}}}}=\tilde{\ensuremath{{\bm{\mathrm{A}}}}}$,
$\ensuremath{{\bm{\mathrm{X}}}}=\ensuremath{{\bm{\mathrm{V}}}}_{\tilde{\ensuremath{{\bm{\mathrm{A}}}}},k}$, $\mat Y=\ensuremath{{\bm{\mathrm{V}}}}_{\tilde{\ensuremath{{\bm{\mathrm{A}}}}},k+}$,
$\mat L=\diag{\tilde{\lambda}_{1},\dots,\tilde{\lambda}_{k}}$, $\ensuremath{{\bm{\mathrm{M}}}}=\diag{\tilde{\lambda}_{k+1},\dots,\tilde{\lambda}_{n}},$
$\mat Z=\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$, $\ensuremath{{\bm{\mathrm{N}}}}=\diag{\lambda_{1},\dots,\lambda_{k}}$,
and $\delta=\lambda_{k}-\tilde{\lambda}_{k+1}$ . It is easy to verify
that the conditions of the $\sin\Theta$ Theorem hold, so
\[
\TNorm{\sin\Theta(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k},\ensuremath{{\bm{\mathrm{V}}}}_{\tilde{\ensuremath{{\bm{\mathrm{A}}}}},k})}\leq\frac{\TNorm{\mat R}}{\lambda_{k}-\tilde{\lambda}_{k+1}}
\]
where $\mat R=\tilde{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}-\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}\ensuremath{{\bm{\mathrm{N}}}}$.
We have $\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}=\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}\ensuremath{{\bm{\mathrm{N}}}}$ so $\TNorm{\mat R}=\TNorm{(\tilde{\ensuremath{{\bm{\mathrm{A}}}}}-\ensuremath{{\bm{\mathrm{A}}}})\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}\leq\TNorm{\tilde{\ensuremath{{\bm{\mathrm{A}}}}}-\ensuremath{{\bm{\mathrm{A}}}}}$.
Combining this inequality with the previous one and noting that $d_{2}(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k},\ensuremath{{\bm{\mathrm{V}}}}_{\tilde{\ensuremath{{\bm{\mathrm{A}}}}},k})=\TNorm{\sin\Theta(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k},\ensuremath{{\bm{\mathrm{V}}}}_{\tilde{\ensuremath{{\bm{\mathrm{A}}}}},k})}$
completes the proof.
\end{proof}
Under the conditions of Theorem~\ref{thm:sin-theta-1}, since $\ensuremath{{\bm{\mathrm{A}}}}$
and $\tilde{\ensuremath{{\bm{\mathrm{A}}}}}$ are symmetric matrices, Weyl's inequality implies
that
\[
d_{2}(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k},\ensuremath{{\bm{\mathrm{V}}}}_{\tilde{\ensuremath{{\bm{\mathrm{A}}}}},k})\leq\frac{\TNorm{\ensuremath{{\bm{\mathrm{A}}}}-\tilde{\ensuremath{{\bm{\mathrm{A}}}}}}}{\lambda_{k}-\lambda_{k+1}-\TNorm{\ensuremath{{\bm{\mathrm{A}}}}-\tilde{\ensuremath{{\bm{\mathrm{A}}}}}}}\,
\]
as long as $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}-\tilde{\ensuremath{{\bm{\mathrm{A}}}}}}<\lambda_{k}-\lambda_{k+1}$.
Thus, if $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}-\tilde{\ensuremath{{\bm{\mathrm{A}}}}}}\ll\lambda_{k}-\lambda_{k+1}$
then we can compute an approximation to the $k$-dimensional dominant
subspace of $\ensuremath{{\bm{\mathrm{A}}}}$ by computing the $k$-dimensional dominant subspace
of $\tilde{\ensuremath{{\bm{\mathrm{A}}}}}$.

\section{\label{sec:dim-reduce-pcr}PCR with Dimensionality Reduction}

Our goal is to design algorithms which compute an approximate solution
to the PCR or PCP problem. Our strategy for designing such algorithms
is to reduce the dimensions of $\ensuremath{{\bm{\mathrm{A}}}}$ prior to computing the PCR/PCP
solution. Specifically, let $\mat R\in\mathbb{R}^{d\times t}$ be some matrix
where $t\leq d$, and define
\begin{equation}
{\mathbf{x}}_{\mat R,k}:=\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}(\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}{\mathbf{b}}\,.\label{eq:apcr}
\end{equation}
The rationale in Eq.~(\ref{eq:apcr}) is as follows. First, $\ensuremath{{\bm{\mathrm{A}}}}$
is compressed by computing $\ensuremath{{\bm{\mathrm{A}}}}\mat R$ (this is the dimensionality
reduction step). Then we compute the rank $k$ PCR solution of $\ensuremath{{\bm{\mathrm{A}}}}\mat R$
and ${\mathbf{b}}$; this is $(\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}{\mathbf{b}}$. Finally,
the solution is projected back to the original space by multiplying
by $\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}$. Obviously, given $\mat R$ we can
compute ${\mathbf{x}}_{\mat R,k}$ in $O(ndt)$ (and even faster, if $\ensuremath{{\bm{\mathrm{A}}}}$
is sparse), so if $t\ll\min(n,d)$ there is a potential for significant
gain in terms of computational complexity\emph{ }provided it is possible
to compute $\mat R$ efficiently as well. Furthermore, if we design
$\mat R$ to have some special structure that allows us to compute
$\ensuremath{{\bm{\mathrm{A}}}}\mat R$ in $O(nt^{2})$ time, the overall complexity would reduce
to $O(nt^{2})$.

Of course, ${\mathbf{x}}_{\mat R,k}$ is not the PCR solution ${\mathbf{x}}_{k}$ (unless
$\mat R=\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$). This suggests the following mathematical
question (which, in turn, leads to an algorithmic question): under
which conditions on $\mat R$ is ${\mathbf{x}}_{\mat R,k}$ a good approximation
to the PCR solution ${\mathbf{x}}_{k}$? In this section, we derive general
conditions on $\mat R$ that ensure deterministically that ${\mathbf{x}}_{\mat R,k}$
is in some sense (which we formalize later in this section) a good
approximation of ${\mathbf{x}}_{k}$. The results in this section are non algorithmic
and independent of the method in which $\mat R$ is computed. In the
next section we address the algorithmic question: how can we compute
such $\mat R$ matrices efficiently?

We approach the mathematical question from two different perspectives:
an optimization perspective and a statistical perspective. In the
optimization perspective, we consider PCR/PCP as an optimization problem
(Eq.~(\ref{eq:pcr-opt-problem})), and ask whether the value of the
objective function of ${\mathbf{x}}_{\mat R,k}$ is close to optimal value of
the objective function, while upholding the constraints approximately
(see Definition~\ref{def:apcr-apcp}). In the statistical perspective,
we treat ${\mathbf{x}}_{k}$ and ${\mathbf{x}}_{\mat R,k}$ as statistical estimators,
and compare their excess risk under a fixed-design model. Interestingly,
the conditions we derive for $\mat R$ are the same for both perspectives.

Before proceeding, we remark that an important special case of (\ref{eq:apcr})
is when $\mat R$ has exactly $k$ columns. In that case, for brevity,
we omit the subscript $k$ from ${\mathbf{x}}_{\mat R,k}$ and notice that
\begin{equation}
{\mathbf{x}}_{\mat R}=\mat R(\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{+}{\mathbf{b}}\,.\label{eq:cls}
\end{equation}
Eq.~(\ref{eq:cls}) is valid even if $\mat R$ has more than $k$
columns and/or the columns are not orthonormal. Thus, an established
technique in the literature, frequently referred to as Compressed
Least Squares (CLS)~\cite{maillard2009compressed,Kaban14,slawski2017compressed,Slawski17,THM17},
is to generate a random $\mat R$ and compute ${\mathbf{x}}_{\mat R}$. To avoid
confusion, we stress the difference between (\ref{eq:apcr}) and (\ref{eq:cls}):
in (\ref{eq:apcr}) we compute a PCR solution on the compressed matrix
$\ensuremath{{\bm{\mathrm{A}}}}\mat R$, while in~(\ref{eq:cls}) ordinary least squares is
used. These two strategies coincide when $\mat R$ has $k$ columns.
In this paper, we focus on Eq.~(\ref{eq:apcr}) and consider Eq.~(\ref{eq:cls})
only when it is a special case of Eq.~(\ref{eq:apcr}) (when $\mat R$
has exactly $k$ columns). For an analysis of CLS from a statistical
perspective, see recent work by Slawski~\cite{Slawski17}.

\subsection{Optimization Perspective}

The PCR solution can be written as the solution of a constrained least
squares problem:
\[
{\mathbf{x}}_{k}=\arg\min_{\begin{array}{c}
\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\textsc{T}}{\mathbf{x}}}=0\\
{\mathbf{x}}\in\range{\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}}
\end{array}}\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}-{\mathbf{b}}}.
\]
In order to analyze a candidate solution $\tilde{{\mathbf{x}}}$ from an optimization
perspective, we need to decide how to treat the constraints. One option
is to require a candidate $\tilde{{\mathbf{x}}}$ to be inside the feasible
set. Indeed, Pilanci and Wainwright recently considered sketching
based methods for constrained least squares regression~\cite{pilanci2015randomized}.
However, there is no evident way to impose $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\mat T}{\mathbf{x}}=0$
without actually computing $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}$, which is as expensive
as computing $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$. Thus, if we require an approximate
solution to be inside the feasible set, we might as well compute the
exact PCR solution. Thus, in our notion of approximate PCR, we relax
the constraints and require only that the approximate solution is
close to meeting the constraint, i.e. we seek a solution for which
$\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\tilde{{\mathbf{x}}}-{\mathbf{b}}}$ is close to $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}-{\mathbf{b}}}$
\emph{and} $\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\textsc{T}}\tilde{{\mathbf{x}}}}$ is small.

Similarly, the if $\ensuremath{{\bm{\mathrm{A}}}}$ has full rank the PCP solution can written
as the solution of a constrained least squares problem:
\[
{\mathbf{b}}_{k}=\arg\min_{\begin{array}{c}
\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\textsc{T}}\tilde{{\mathbf{b}}}}=0\\
\tilde{{\mathbf{b}}}\in\range{\ensuremath{{\bm{\mathrm{A}}}}}
\end{array}}\TNorm{\tilde{{\mathbf{b}}}-{\mathbf{b}}}
\]
Again, our notion of approximate PCP relaxes the constraint.

The discussion above motivates the following definition of approximate
PCR/PCP.
\begin{defn}
[Approximate PCR and PCP] \label{def:apcr-apcp}An estimator $\tilde{{\mathbf{x}}}$
is an \emph{$(\epsilon,\upsilon)$-approximate PCR} of rank $k$ if
\[
\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\tilde{{\mathbf{x}}}-{\mathbf{b}}}=\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}-{\mathbf{b}}}\pm\epsilon\TNorm{{\mathbf{b}}}
\]
and $\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\textsc{T}}\tilde{{\mathbf{x}}}}\leq\upsilon\TNorm{{\mathbf{b}}}$.
An estimator $\tilde{{\mathbf{b}}}$ is an \emph{($\epsilon,\upsilon)$-approximate
PCP} of rank $k$ if
\[
\TNorm{\tilde{{\mathbf{b}}}-{\mathbf{b}}}=\TNorm{{\mathbf{b}}_{k}-{\mathbf{b}}}\pm\epsilon\TNorm{{\mathbf{b}}}
\]
and $\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\textsc{T}}\tilde{{\mathbf{b}}}}\leq\upsilon\TNorm{{\mathbf{b}}}$.
\end{defn}

Before proceeding, a few remarks are in order.
\begin{enumerate}
\item Imposing no constraints on $\tilde{{\mathbf{x}}}$ (or $\tilde{{\mathbf{b}}}$) does not
make sense: we can always form or approximate the ordinary least squares
solution and it will demonstrate a smaller objective value. Indeed,
the main motivation for using PCR to impose some form of regularization,
so it is crucial the definition of approximate PCR/PCP have some form
of regularization built-in.
\item We require only additive error on the objective function, while relative
error bounds are usually viewed as more desirable. For approximate
PCR, requiring relative error bounds is likely unrealistic: since
it is possible that ${\mathbf{b}}=\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}$, any algorithm that provides
a relative error bound must search inside a space that contains $\range{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}$.
This is a strong restriction (and plausibly one that actually requires
computing $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$ ).
\item Approximate PCR implies approximate PCP: if $\tilde{{\mathbf{x}}}$ is an $(\epsilon,\nu)$-approximate
PCR then $\ensuremath{{\bm{\mathrm{A}}}}\tilde{{\mathbf{x}}}$ is an $(\epsilon,\sigma_{k+1}\nu)$-approximate
PCP.
\item Our notion of approximate PCP is somewhat similar to the notion of
approximate PCP proposed recently by Allen-Zhu and Li~\cite{AL17}.
\item Yet another notion of approximate PCR appears in~\cite[Theorem 5]{BM14}.
They too, consider an additive error on objective function, but instead
of considering the distance to the dominant subspace they bound the
distance of the approximate solution to the true solution. We remark
that a bound on $\TNorm{{\mathbf{x}}_{k}-\tilde{{\mathbf{x}}}}$ trivially implies a bound
on $\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\textsc{T}}\tilde{{\mathbf{x}}}}$.
\item Arguably, it would have been preferable to require the approximate
PCR solution $\tilde{{\mathbf{x}}}$ to be such that $\TNorm{{\mathbf{x}}_{k}-\tilde{{\mathbf{x}}}}$
is small (relative to $\TNorm{{\mathbf{x}}_{k}}$). However, be believe that
providing such guarantees with reasonable sketch sizes requires iterations.
In this paper, we focus predominately on algorithms that do not require
iterations (the only exception being the input sparsity algorithm
in subsection~\ref{subsec:input-sparsity}).
\end{enumerate}
We are now ready to state general conditions on $\mat R$ that ensure
deterministically that ${\mathbf{x}}_{\mat R}$ is an approximate PCR, and conditions
on $\mat R$ that ensure deterministically that $\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}$
is an approximate PCP.

\begin{thm}
\label{thm:structural}Suppose that $\mat R\in\mathbb{R}^{d\times s}$ where
$s\ge k$. Assume that $\nu\in(0,1)$.
\begin{enumerate}
\item If $d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\leq\nu$
then $\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}$ is an $(\nu,\nu)$-approximate PCP.
\item \label{enu:structural_k_equal_s}If $s=k$, $\mat R$ has orthonormal
columns (i.e., $\mat R^{\textsc{T}}\mat R=\mat I_{k})$ and $d_{2}\left(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\right)\leq\nu(1+\nu^{2})^{-1/2}$
then ${\mathbf{x}}_{\mat R}$ is an $\left(\frac{\sigma_{k+1}}{\sigma_{k}}\nu,\frac{\nu}{\left(\sqrt{1-\nu^{2}}-\nu\right)\sigma_{k}}\right)$-approximate
PCR.
\end{enumerate}
\end{thm}

Before proving this theorem, we state a theorem which is a corollary
of a more general result proved recently by Drineas et al.~\cite{drineas2016structural},
and then proceed to proving a couple of auxiliary lemmas.
\begin{thm}
[Corollary of Theorem 2.1 in \cite{drineas2016structural}]\label{thm:krylov_srtuctural}
Let $\ensuremath{{\bm{\mathrm{A}}}}$ be an $m\times n$ matrix with singular value decomposition
$\ensuremath{{\bm{\mathrm{A}}}}=\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}^{\mat T}$ .
Let $k\ge0$ and let $\mat R\in\mathbb{R}^{d\times k}$ be any matrix such
that $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\mat T}\mat R$ has full rank. Then,
\end{thm}

\[
\TNorm{\sin\Theta(\ensuremath{{\bm{\mathrm{A}}}}\mat R,\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})}\leq\TNorm{\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}},k+}}\cdot\TNorm{\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{-1}}\cdot\TNorm{\tan\Theta(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})}
\]
\begin{lem}
\label{lem:11}Assume $\rank{\ensuremath{{\bm{\mathrm{A}}}}}\ge k$. If $\mat R\in\mathbb{R}^{d\times k}$
has orthonormal columns and $d_{2}\left(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\right)\leq\nu$
then the following bounds hold:

\begin{equation}
\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R}\leq\nu\label{eq:orho-to-minor}
\end{equation}
\begin{equation}
\sigma_{\min}\left(\ensuremath{{\bm{\mathrm{A}}}}\mat R\right)\ge\sigma_{k}\left(\sqrt{1-\nu^{2}}-\nu\right)\label{eq:singular_values_bound}
\end{equation}
Furthermore, if $\nu<1$ then $\rank{\ensuremath{{\bm{\mathrm{A}}}}\mat R}=k$.
\end{lem}

\begin{proof}
Since both $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$ and $\mat R$ have orthonormal columns,
$d_{2}\left(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\right)\leq\nu$ implies that the
square of the singular values of $\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}\mat R$ lie
inside the interval $[1-\nu^{2},1]$. The eigenvalues of $\mat R^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}\mat R$
are exactly the square of the singular values of $\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}\mat R$,
so the eigenvalues of $\mat I_{k}-\mat R^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}\mat R$
lie in $[0,\nu^{2}]$. Let $\mat Z$ be any matrix with orthonormal
columns that completes $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}$ to a basis (i.e. $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}^{\textsc{T}}+\mat Z\matZ^{\textsc{T}}=\mat I_{d}$)
and is orthogonal to $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}$ (i.e., $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}^{\textsc{T}}\mat Z=0$).
Note that $\mat Z$ can be an empty matrix if $d\leq n$. Denote $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k\perp}=\left[\begin{array}{cc}
\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+} & \mat Z\end{array}\right]$. We have
\begin{eqnarray*}
\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k\perp}^{\textsc{T}}\mat R} & = & \sqrt{\TNorm{\mat R^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k\perp}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k\perp}^{\textsc{T}}\mat R}}\\
 & = & \sqrt{\TNorm{\mat I_{k}-\mat R^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}\mat R}}\\
 & \leq & \nu
\end{eqnarray*}
where we used the fact that $\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}+\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k\perp}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k\perp}^{\textsc{T}}=\mat I_{d}$.
We now note that $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\textsc{T}}\mat R$ is a submatrix of $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k\perp}^{\textsc{T}}$
so $\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R}\leq\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k\perp}^{\textsc{T}}\mat R}\leq\nu$.
This establishes the first part of the theorem.

As for the second part, recall the following identities: 1) for any
matrix $\ensuremath{{\bm{\mathrm{X}}}}$ and $\mat Y$ of the same size: $\sigma_{\min}\left(\ensuremath{{\bm{\mathrm{X}}}}\pm\mat Y\right)\ge\sigma_{\min}\left(\ensuremath{{\bm{\mathrm{X}}}}\right)-\sigma_{\max}\left(\mat Y\right)$
\cite[Theorem 3.3.19]{horn1990matrix}, 2) if the number of rows in
$\ensuremath{{\bm{\mathrm{X}}}}$ and $\mat Y$ is at least as large as the number of columns,
and $\ensuremath{{\bm{\mathrm{X}}}}\mat Y$ is defined, then $\sigma_{\min}\left(\ensuremath{{\bm{\mathrm{X}}}}\mat Y\right)\ge\sigma_{\min}\left(\ensuremath{{\bm{\mathrm{X}}}}\right)\sigma_{\min}\left(\mat Y\right)$.
We have

\begin{eqnarray*}
\sigma_{\min}\left(\ensuremath{{\bm{\mathrm{A}}}}\mat R\right) & = & \sigma_{\min}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}\mat R+\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R)\\
 & \geq & \sigma_{\min}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}\mat R)-\sigma_{\max}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R)\\
 & \geq & \sigma_{\min}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k})\sigma_{\min}(\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}\mat R)-\sigma_{\max}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+})\sigma_{\max}(\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R)\\
 & = & \sigma_{k}\sigma_{\min}(\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}\mat R)-\sigma_{k+1}\sigma_{\max}(\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R)\\
 & \geq & \sigma_{k}\sqrt{1-\nu^{2}}-\sigma_{k+1}\nu\\
 & \geq & \sigma_{k}(\sqrt{1-\nu^{2}}-\nu)
\end{eqnarray*}
where the first equality follows from the fact that $\ensuremath{{\bm{\mathrm{A}}}}(\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}^{\textsc{T}}+\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}})=\ensuremath{{\bm{\mathrm{A}}}}$.
When $\nu<1$ we have $\sigma_{\min}(\ensuremath{{\bm{\mathrm{A}}}}\mat R)>0$, so indeed the
rank of $\ensuremath{{\bm{\mathrm{A}}}}\mat R$ is $k$.
\end{proof}
\begin{lem}
\label{lem:14}Assume $\rank{\ensuremath{{\bm{\mathrm{A}}}}}\ge k$. Suppose that $\mat R\in\mathbb{R}^{d\times k}$
has orthonormal columns, and that $d_{2}\left(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\right)\leq\nu(1+\nu^{2})^{-1/2}<1$.
We have
\[
d_{2}(\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R},\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})\leq\frac{\sigma_{k+1}}{\sigma_{k}}\nu\,.
\]
\end{lem}

\begin{proof}
Since $\rank{\ensuremath{{\bm{\mathrm{A}}}}}\ge k$ and $\nu(1+\nu^{2})^{-1/2}<1$, according
to Lemma \ref{lem:11} the matrix $\ensuremath{{\bm{\mathrm{A}}}}\mat R$ has full rank. According
to Theorem~ \ref{thm:projection_equality} and the fact that $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R}$
and $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}$ are orthogonal projections we have
\begin{equation}
d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R},\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}\right)=d_{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\mat R,\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}\right)=\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}}=\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R}^{\perp}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}}\label{eq:projection_equality_k0}
\end{equation}
Combining Theorem \ref{thm:krylov_srtuctural} and Eq. (\ref{eq:projection_equality_k0}),
we bound:

\begin{eqnarray*}
d_{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\mat R,\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}\right) & = & \TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R}^{\perp}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}}\\
 & = & \TNorm{(\mat I-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R})\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}\\
 & = & \TNorm{\sin\Theta(\ensuremath{{\bm{\mathrm{A}}}}\mat R,\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})}\\
 & \leq & \TNorm{\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}.k+}}\cdot\TNorm{\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}.k}^{-1}}\cdot\TNorm{\tan\Theta(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})}\\
 & = & \frac{\sigma_{k+1}}{\sigma_{k}}\cdot\TNorm{\tan\Theta(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})}\\
 & \leq & \frac{\sigma_{k+1}}{\sigma_{k}}\nu
\end{eqnarray*}
where the last inequality follows from the fact that $\Theta(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})$
is a diagonal matrix whose diagonal values are the inverse cosine
of the singular values of $\mat R^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$, and these,
in turn, are all larger than $\sqrt{1-\nu^{2}(1+\nu^{2})^{-1}}.$
\end{proof}
We are now ready to prove Theorem~\ref{thm:structural}.
\begin{proof}
[Proof of Theorem~\ref{thm:structural}]We need to show both the additive
error bounds on the objective function, and the error bound on the
constraints. We start with the additive error bounds on the objective
function, both for PCP (first part of the theorem) and for PCR (second
part of the theorem). We have
\[
\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}=\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\left(\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\right)^{+}{\mathbf{b}}=\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k}}{\mathbf{b}}
\]
and
\[
\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}=\mat A\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})^{+}{\mathbf{b}}=\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}}{\mathbf{b}}\,.
\]
Thus,
\begin{eqnarray*}
\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}-{\mathbf{b}}} & = & \TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}-{\mathbf{b}}+\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}-\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}}\\
 & = & \TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}-{\mathbf{b}}}\pm\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}-\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}}\\
 & = & \TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}-{\mathbf{b}}}\pm\TNorm{(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}){\mathbf{b}}}\\
 & = & \TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}-{\mathbf{b}}}\pm d_{2}(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k})\cdot\TNorm{{\mathbf{b}}}
\end{eqnarray*}
In the first part of the theorem, we have $d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\leq\nu$,
while the second part of the theorem we have $\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}=\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R}$
(since $\ensuremath{{\bm{\mathrm{A}}}}\mat R$ has $k$ columns) and Lemma~\ref{lem:14} ensures
that $d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\leq\nu\sigma_{k+1}/\sigma_{k}$.
Either way, the additive error bounds of the theorem are met.

We now bound the infeasibility of the approximate solution for the
PCP guarantee (first part of the theorem):
\begin{eqnarray*}
\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}} & = & \TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}-\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}+\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}}\\
 & \leq & \TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\left(\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}-\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}\right)}+\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}}\\
 & \leq & \TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}-\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}}\\
 & \leq & d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\TNorm{{\mathbf{b}}}\\
 & \leq & \nu\TNorm{{\mathbf{b}}}
\end{eqnarray*}
where we used the fact that $\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}\in\range{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}$
so $\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}=0$.

We now bound the infeasibility of the approximate solution for the
PCR guarantee (second part of the theorem):

\begin{eqnarray*}
\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}{\mathbf{x}}_{\mat R}} & = & \TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R(\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{+}{\mathbf{b}}}\\
 & \leq & \TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R}\cdot\TNorm{(\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{+}}\cdot\TNorm{{\mathbf{b}}}\\
 & \leq & \frac{\nu}{\left(\sqrt{1-\nu^{2}}-\nu\right)\sigma_{k}}\TNorm{{\mathbf{b}}}
\end{eqnarray*}
where we used Lemma~\ref{lem:11} to bound $\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R}$
and $\TNorm{(\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{+}}$.
\end{proof}

\subsection{\label{subsec:statistical}Statistical Perspective}

We now consider ${\mathbf{x}}_{\mat R,k}$ from a statistical perspective. We
use a similar framework to the one used in the literature to analyze
CLS~\cite{slawski2017compressed,Slawski17,THM17}. That is, we consider
a fixed design setting in which the rows of $\ensuremath{{\bm{\mathrm{A}}}}$, ${\mathbf{a}}_{1},\dots,{\mathbf{a}}_{n}\in\mathbb{R}^{d}$,
are considered as fixed, and ${\mathbf{b}}$'s entries, $b_{1},\dots,b_{n}\in\mathbb{R}$,
are
\[
b_{i}=f_{i}+\xi_{i}
\]
where $f_{1},\dots,f_{n}$ are fixed values and the noise terms $\xi_{1},\dots,\xi_{n}$
are assumed to be independent random values with zero mean and $\sigma^{2}$
variance. We denote by ${\mathbf{f}}\in\mathbb{R}^{n}$ the vector whose $i$th entry
is $f_{i}$ . The goal is to recover ${\mathbf{f}}$ from ${\mathbf{b}}$ (i.e., de-noise
${\mathbf{b}}$).

The optimal predictor $\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}^{\star}$ of ${\mathbf{f}}$ given $\ensuremath{{\bm{\mathrm{A}}}}$ is
a minimizer of
\[
\min_{{\mathbf{x}}\in\mathbb{R}^{d}}\Expect{\TNormS{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}-{\mathbf{b}}}/n}
\]
where here, and in subsequent expressions, the expectation is with
respect to the noise $\xi$ (if there are multiple minimizers, ${\mathbf{x}}^{\star}$
is the minimizer with minimum norm). It is easy to verify that $\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}^{\star}=\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}{\mathbf{f}}$.

Given an estimator $\theta=\theta(\ensuremath{{\bm{\mathrm{A}}}},{\mathbf{b}})$ of ${\mathbf{x}}^{\star}$ (which
we assume is a random variable since ${\mathbf{b}}$ is a random varaible),
its \emph{excess risk} is define as
\[
{\cal E}(\theta)\coloneqq\Expect{\TNormS{\ensuremath{{\bm{\mathrm{A}}}}\theta-\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}^{\star}}/n}\,.
\]
The \emph{ordinary least square estimator }(OLS) $\hat{{\mathbf{x}}}$ is simply
a solution to $\min_{{\mathbf{x}}\in\mathbb{R}^{d}}\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}-{\mathbf{b}}}$: $\hat{{\mathbf{x}}}\coloneqq\ensuremath{{\bm{\mathrm{A}}}}^{+}{\mathbf{b}}$.
Simple calculations show that
\[
{\cal E}(\hat{{\mathbf{x}}})=\sigma^{2}\rank{\ensuremath{{\bm{\mathrm{A}}}}}/n\,.
\]
Thus, if the rank of $\ensuremath{{\bm{\mathrm{A}}}}$ is large, which is usually the case
when $d\gg n$, then the excess risk might be large (and it does not
asymptotically converge to $0$ if $\rank{\ensuremath{{\bm{\mathrm{A}}}}}=\Omega(n)$). This
motivates the use of regularization (e.g., PCR). Indeed, the excess
risk of the PCR estimator ${\mathbf{x}}_{k}$ can be bounded~\cite{Slawski17}:
\begin{equation}
{\cal E}({\mathbf{x}}_{k})\leq\frac{\InfNormS{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}^{\textsc{T}}{\mathbf{x}}^{\star}}\cdot\sum_{i=k+1}^{\min(n,d)}\sigma_{i}^{2}}{n}+\frac{\sigma^{2}k}{n}\,.\label{eq:pcr-old-risk-bound}
\end{equation}
In many scenarios, ${\mathbf{x}}_{k}$ has a significantly reduced excess risk
in comparison to the excess risk of $\hat{{\mathbf{x}}}$ (see~\cite{Slawski17}
for a discussion). This motivates the use of PCR when $d$ is large.

In this section, we analyze the excess risk of ${\mathbf{x}}_{\mat R,k}$ based
on properties of $\mat R$. The bounds are based on the following identity~\cite{Slawski17}\footnote{However, no proof of (\ref{eq:bias-var}) appears in~\cite{Slawski17},
so for completeness we include a proof in the appendix. }: for any $\ensuremath{{\bm{\mathrm{M}}}}$ of appropriate size
\begin{equation}
{\cal E}({\mathbf{x}}_{\ensuremath{{\bm{\mathrm{M}}}}})={\cal E}(\ensuremath{{\bm{\mathrm{M}}}}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{M}}}})^{\textsc{+}}{\mathbf{b}})=\underset{{\cal B}({\mathbf{x}}_{\ensuremath{{\bm{\mathrm{M}}}}})}{\underbrace{\frac{1}{n}\TNormS{(\mat I-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{M}}}}})\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}^{\star}}}}+\underset{{\cal V}({\mathbf{x}}_{\ensuremath{{\bm{\mathrm{M}}}}})}{\underbrace{\sigma^{2}\frac{\rank{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{M}}}}}}{n}}}\,.\label{eq:bias-var}
\end{equation}
In the above, ${\cal B}({\mathbf{x}}_{\ensuremath{{\bm{\mathrm{M}}}}})$ can be viewed as a bias term,
and ${\cal V}({\mathbf{x}}_{\ensuremath{{\bm{\mathrm{M}}}}})$ can be viewed as a variance term. Eq.~(\ref{eq:pcr-old-risk-bound})
is obtained by bounding the bias term ${\cal B}({\mathbf{x}}_{k})$, although
our results lead to a bound on ${\cal E}({\mathbf{x}}_{k})$ that is tighter
in some cases (Corollary~\ref{cor:new-pcr-risk}). An immediate corollary
of~(\ref{eq:bias-var}) is the following bound for ${\mathbf{x}}_{\mat R,k}$:
\begin{equation}
{\cal E}({\mathbf{x}}_{\mat R,k})=\frac{1}{n}\TNormS{(\mat I-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}})\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}^{\star}}+\frac{\sigma^{2}k}{n}\,.\label{eq:pcp_risk}
\end{equation}

The following results addresses the case where $\mat R$ has $k$ orthonormal
columns. The conditions are the same as the first part of Theorem~\ref{thm:structural}
(optimization perspective analysis).
\begin{thm}
\label{thm:stat-structural}Assume that $\rank{\ensuremath{{\bm{\mathrm{A}}}}}\geq k$. Suppose
that $\mat R\in\mathbb{R}^{d\times k}$ has orthonormal columns, and that $d_{2}(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})\leq\nu(1+\nu^{2})^{-1/2}<1$.
Then,
\[
{\cal E}({\mathbf{x}}_{\mat R})\leq\frac{\left(1+\nu\right)\cdot\TNormS{{\mathbf{x}}^{\star}}\cdot\sigma_{k+1}^{2}}{n}+\frac{\sigma^{2}k}{n}
\]
\end{thm}

For the proof, we need the following theorem due to Halko et al.~\cite{halko2011finding}.
\begin{thm}
[Theorem 9.1 in \cite{halko2011finding}]\label{thm:hmt_structural}
Let $\ensuremath{{\bm{\mathrm{A}}}}$ be an $m\times n$ matrix with singular value decomposition
$\ensuremath{{\bm{\mathrm{A}}}}=\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}^{\mat T}$ .Let
$k\ge0$. Let $\mat R$ be any matrix such that $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\mat T}\mat R$
has full row rank. Then we have
\[
\TNormS{(\mat I_{m}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R})\ensuremath{{\bm{\mathrm{A}}}}}\leq\TNormS{\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}},k+}}+\TNormS{\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}},k+}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\textsc{T}}\mat R\left(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\textsc{T}}\mat R\right)^{+}}
\]
\end{thm}

\begin{proof}
[Proof of Theorem~\ref{thm:stat-structural}]The condition that $d_{2}(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})\leq\nu(1+\nu^{2})^{-1/2}<1$
ensures that $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\textsc{T}}\mat R$ has full rank, and that
$\TNormS{\tan\Theta(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})}\leq\nu$ (since $\Theta(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})$
is a diagonal matrix whose diagonal values are the inverse cosine
of the singular values of $\mat R^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}$, and these, in turn,
are all larger than $\sqrt{1-\nu^{2}(1+\nu^{2})^{-1}}$). Thus we
have,

\begin{eqnarray*}
{\cal B}({\mathbf{x}}_{\mat R}) & = & \frac{1}{n}\TNormS{(\mat I-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R})\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}^{\star}}\\
 & \leq & \frac{1}{n}\TNormS{{\mathbf{x}}^{\star}}\cdot\left(\TNormS{\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}},k+}}+\TNormS{\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}},k+}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\mat T}\mat R\left(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\mat T}\mat R\right)^{+}}\right)\\
 & \leq & \frac{1}{n}\TNormS{{\mathbf{x}}^{\star}}\left(\sigma_{k+1}^{2}+\sigma_{k+1}^{2}\TNormS{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k+}^{\mat T}\mat R\left(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\mat T}\mat R\right)^{+}}\right)\\
 & = & \frac{1}{n}\TNormS{{\mathbf{x}}^{\star}}\left(\sigma_{k+1}^{2}+\sigma_{k+1}^{2}\TNormS{\tan\Theta\left(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}\right)}\right)\\
 & \leq & \frac{\left(1+\nu\right)\cdot\TNormS{{\mathbf{x}}^{\star}}\cdot\sigma_{k+1}^{2}}{n}
\end{eqnarray*}
where in the first inequality we used Theorem~\ref{thm:hmt_structural}
and for the second equality we used Lemma~\ref{lem:tan-to-spectral}.
The result now follows from the fact that $\rank{\ensuremath{{\bm{\mathrm{A}}}}\mat R}\leq k$.
\end{proof}
\begin{cor}
\label{cor:new-pcr-risk}For the PCR solution ${\mathbf{x}}_{k}$ we have
\[
{\cal E}({\mathbf{x}}_{k})\leq\frac{\TNormS{{\mathbf{x}}^{\star}}\cdot\sigma_{k+1}^{2}}{n}+\frac{\sigma^{2}k}{n}\,.
\]
\end{cor}

Next, we consider the general case where $\mat R$ does not necessarily
have orthonormal columns, and potentially has more than $k$ columns.
The conditions are the same as the second part of Theorem~\ref{thm:structural}
(optimization perspective).
\begin{thm}
\label{thm:struct-stat-pcp}Suppose that $\mat R\in\mathbb{R}^{d\times s}$
where $s\ge k$. Assume that $\rank{\ensuremath{{\bm{\mathrm{A}}}}\mat R}\geq k$. If $d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\leq\nu<1$
then,
\[
{\cal E}({\mathbf{x}}_{\mat R,k})\leq{\cal E}({\mathbf{x}}_{k})+\frac{(2\nu+\nu^{2})\TNormS{{\mathbf{f}}}}{n}\,.
\]
\end{thm}

\begin{proof}
Since $\ensuremath{{\bm{\mathrm{A}}}}\mat R$ has rank at least $k$, we have $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}=\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}$.
From (\ref{eq:pcp_risk}), the fact that $\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}^{\star}=\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}{\mathbf{f}}$,
and $\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}=\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}$
(since the range of $\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}$ is contained
in the range of $\ensuremath{{\bm{\mathrm{A}}}}$) we have

\begin{eqnarray*}
{\cal B}({\mathbf{x}}_{\mat R,k}) & = & \frac{1}{n}\TNormS{(\mat I-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}})\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}^{\star}}\\
 & = & \frac{1}{n}\TNormS{(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}){\mathbf{f}}}\\
 & = & \frac{1}{n}\TNormS{(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}+\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}){\mathbf{f}}}\\
 & = & \frac{1}{n}\left(\TNormS{(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}){\mathbf{f}}}+\TNormS{(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}){\mathbf{f}}}+2(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}{\mathbf{f}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}{\mathbf{f}})^{\textsc{T}}(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}{\mathbf{f}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}{\mathbf{f}})\right)\\
 & \leq & {\cal B}({\mathbf{x}}_{k})+d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)^{2}\frac{\TNormS{{\mathbf{f}}}}{n}+\frac{2}{n}\left|{\mathbf{f}}^{\textsc{T}}(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}})^{\textsc{T}}(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}){\mathbf{f}}\right|
\end{eqnarray*}
For the cross-terms, we bound

\begin{eqnarray*}
\left|{\mathbf{f}}^{\textsc{T}}(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}})^{\textsc{T}}(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}){\mathbf{f}}\right| & = & \left|{\mathbf{f}}^{\textsc{T}}\left(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}+\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}\right){\mathbf{f}}\right|\\
 & = & \left|{\mathbf{f}}^{\textsc{T}}\left(\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}+\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}\right){\mathbf{f}}\right|\\
 & = & {\mathbf{f}}^{\textsc{T}}\left(\mat I-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}\right)\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}{\mathbf{f}}\\
 & = & {\mathbf{f}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}^{\perp}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}{\mathbf{f}}\\
 & \leq & \TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}^{\perp}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}}\cdot\TNormS{{\mathbf{f}}}
\end{eqnarray*}
Since both $\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k}$ and $\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}$ are full
rank,\textcolor{red}{{} }we have (Theorem~\ref{thm:projection_equality})
\[
\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}^{\perp}\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}}=\TNorm{\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}}-\ensuremath{{\bm{\mathrm{P}}}}_{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}}=d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)
\]
Thus, we find that
\[
{\cal B}({\mathbf{x}}_{\mat R,k})\leq{\cal B}({\mathbf{x}}_{k})+(2\nu+\nu^{2})\frac{\TNormS{{\mathbf{f}}}}{n}\,.
\]
We reach the bound in the theorem statement by adding the variance
${\cal V}({\mathbf{x}}_{\mat R,k})$, which is equal to the variance of ${\mathbf{x}}_{k}$
because the ranks are equal.
\end{proof}

\paragraph{Discussion. }

Theorem~\ref{thm:stat-structural} shows that if $\mat R$ is a good
approximation to $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$, then there is a small relative
increase to the bias term, while the variance term does not change.
Since we are mainly interested in keeping the asymptotic behavior
of the excess risk (as $n$ goes to infinity), a fixed $\nu$ of modest
value suffices. However, for this result to hold, $\mat R$ has to
have exactly $k$ columns and those columns should be orthonormal.
Without these restrictions, we need to resort to Theorem~\ref{thm:struct-stat-pcp}.
In that theorem, we get (if the conditions are met) only an additive
increase in the bias term. Thus if, for example, $\TNormS{{\mathbf{f}}}/n\to c$
as $n\to\infty$ for some constant $c$, then $\nu$ should tend to
$0$ as $n$ goes to infinity, but a constant value should suffice
if $n$ is fixed.

\section{\label{sec:alg}Sketched PCR and PCP}

In the previous section, we considered general conditions on $\mat R$
which ensure that ${\mathbf{x}}_{\mat R,k}$ is an approximate solution to the
PCR/PCP problem. In this section, we propose algorithms to generate
$\mat R$ for which these conditions hold. The main technique we employ
is \emph{matrix sketching}. The idea is to first multiply the data
matrix $\ensuremath{{\bm{\mathrm{A}}}}$ by some random transformation (e.g., a random projection),
and extract an approximate subspace from the compressed matrix.

\subsection{Dimensionality Reduction using Sketching }

The compression (multiplication by a random matrix) alluded in the
previous paragraph can be applied either from the left side, or the
right side, or both. In left sketching, which is more appropriate
if the input matrix has many rows and a modest amount of columns,
we propose to use $\mat R=\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}},k}$ where $\mat S$ is
some sketching matrix (we discuss a couple of options shortly). In
right sketching, which is more appropriate if the input matrix has
many columns and a modest amount of rows, we propose to use $\mat R=\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}$
where $\ensuremath{{\bm{\mathrm{G}}}}$ is some sketching matrix. Two sided sketching, $\mat R=\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k}$
, is aimed for the case that the number of columns and the number
of rows are large.

The sketching matrices, $\mat S$ and $\ensuremath{{\bm{\mathrm{G}}}}$, are randomized dimensionality
reduction transformations. Quite a few sketching transforms have been
proposed in the literature in recent years. For concreteness, we consider
two specific cases, though our results hold for other sketching transformations
as well (albeit some modifications in the bounds might be necessary).
The first, which we refer to as 'subgaussian map', is a random matrix
in which every entry of the matrix is sampled i.i.d from some subgaussian
distribution (e.g. $N(0,1)$) and the matrix is appropriately scaled
(however, scaling is not necessary in our case). The second transform
is a sparse embedding matrix, in which each column is sampled uniformly
and independently from the set of scaled identity vectors and multiplied
by a random sign. We refer to such a matrix as a \noun{CountSketch}
matrix~\cite{charikar2002finding,woodruff2014sketching}.

Both transformations described above, and a few other, have, provided
enough rows are used, with high probability the following property
which we refer to as \emph{approximate Gram property}.
\begin{defn}
Let $\ensuremath{{\bm{\mathrm{X}}}}\in\mathbb{R}^{m\times n}$ be a fixed matrix. For $\epsilon,\delta\in(0,1/2)$,
a distribution ${\cal D}$ on matrices with $m$ columns has the \emph{$(\epsilon,\delta)$-approximate
Gram matrix} property for $\ensuremath{{\bm{\mathrm{X}}}}$ if
\[
\Pr_{\mat S\sim{\cal D}}\left(\TNorm{\ensuremath{{\bm{\mathrm{X}}}}^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{X}}}}-\ensuremath{{\bm{\mathrm{X}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{X}}}}}\geq\epsilon\TNormS{\ensuremath{{\bm{\mathrm{X}}}}}\right)\leq\delta\,.
\]
\end{defn}

Recent results by Cohen et al.~\cite{cohen2015optimal}\footnote{Theorem 1 in \cite{cohen2015optimal} with $k=\sr{\ensuremath{{\bm{\mathrm{X}}}}}.$}
show that when $\mat S$ has independent subgaussian entries, then
as long as the number of rows in $\mat S$ is $\Omega((\sr{\ensuremath{{\bm{\mathrm{X}}}}}+\log(1/\delta))/\epsilon^{2})$
then we have $(\epsilon,\delta)$-approximate Gram property for $\ensuremath{{\bm{\mathrm{X}}}}$.
If $\mat S$ is a \noun{CountSketch} matrix, then as long as the number
of rows in $\mat S$ is $\Omega(\sr{\ensuremath{{\bm{\mathrm{X}}}}}^{2}/(\epsilon^{2}\delta))$
then we have $(\epsilon,\delta)$-approximate Gram property for $\ensuremath{{\bm{\mathrm{X}}}}$~\cite{cohen2015optimal}.

We first describe our results for the various modes of sketching,
and then discuss algorithmic issues and computational complexity.
\begin{thm}
[Left Sketching] Let $\nu,\delta\in(0,1/2)$ and denote
\[
\epsilon=\frac{\nu(1+\nu^{2})^{-1/2}}{1+\nu(1+\nu^{2})^{-1/2}}\cdot\gap{\ensuremath{{\bm{\mathrm{A}}}}}k.
\]
Suppose that $\mat S$ is sampled from a distribution that provides
a $(\epsilon,\delta)$-approximate Gram matrix for $\ensuremath{{\bm{\mathrm{A}}}}$. Then
for $\mat R=\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}},k}$, with probability $1-\delta$, the
approximate solution ${\mathbf{x}}_{\mat R}$ is a $\left(\frac{\sigma_{k+1}}{\sigma_{k}}\nu,\frac{\nu}{\left(\sqrt{1-\nu^{2}}-\nu\right)\sigma_{k}}\right)$-approximate
PCR and
\[
{\cal E}({\mathbf{x}}_{\mat R})\leq\frac{\left(1+\nu\right)\cdot\TNormS{{\mathbf{x}}^{\star}}\cdot\sigma_{k+1}^{2}}{n}+\frac{\sigma^{2}k}{n}
\]

Thus if, for example, $\mat S$ is a \noun{CountSketch} matrix, then
the conditions are met when the number of rows in $\mat S$ is
\[
\Omega\left(\frac{\sr{\ensuremath{{\bm{\mathrm{A}}}}}^{2}}{\gap{\ensuremath{{\bm{\mathrm{A}}}}}k^{2}\nu^{2}\delta}\right)
\]
 rows. In another example, if $\mat S$ is a subgaussian map, then
the conditions are met when the number of rows in $\mat S$ is
\[
\Omega\left(\frac{\sr{\ensuremath{{\bm{\mathrm{A}}}}}+\log(1/\delta)}{\gap{\ensuremath{{\bm{\mathrm{A}}}}}k^{2}\nu^{2}}\right)\,.
\]
\end{thm}

\begin{proof}
Due to Theorems~\ref{thm:structural} and~\ref{thm:stat-structural},
it suffices to show that that $d_{2}\left(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\right)\leq\nu(1+\nu^{2})^{-1/2}$.
Under the conditions of the theorem, with probability of at least
$1-\delta$ we have $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}-\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}}\leq\epsilon\TNormS{\ensuremath{{\bm{\mathrm{A}}}}}$.
If that is indeed the case, $\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}$ has
rank at least $k$ since $\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}$ and $\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}$
are symmetric matrices and we know that $\sigma_{i}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}\right)=\sigma_{i}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}\right)\pm\TNorm{\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}-\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}}$
(Weyl's Theorem and the fact that $\sigma_{k}^{2}>\epsilon\sigma_{1}^{2}$).
Furthermore, since $\nu>0$ we have $\epsilon<\gap{\ensuremath{{\bm{\mathrm{A}}}}}k$, and
Theorem~\ref{thm:sin-theta-1} implies that
\begin{eqnarray*}
d_{2}(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}) & \leq & \frac{\TNorm{\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}-\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}}}{(\sigma_{k}^{2}-\sigma_{k+1}^{2})-\TNorm{\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}-\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}}}\\
 & \leq & \frac{\epsilon}{\gap{\ensuremath{{\bm{\mathrm{A}}}}}k-\epsilon}\\
 & \leq & \nu(1+\nu^{2})^{-1/2}\,.
\end{eqnarray*}
Thus, we have shown that with probability $1-\delta$ we have $d_{2}\left(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\right)\leq\nu(1+\nu^{2})^{-1/2}$,
as required.
\end{proof}
\begin{thm}
[Right Sketching]\label{thm:right_sketching}Let $\nu,\delta\in(0,1/2)$
and denote
\[
\epsilon=\frac{\nu}{1+\nu}\cdot\gap{\ensuremath{{\bm{\mathrm{A}}}}}k.
\]
Suppose that $\ensuremath{{\bm{\mathrm{G}}}}$ is sampled from a distribution that provides
a $(\epsilon,\delta)$-approximate Gram matrix for $\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}$.
Then for $\mat R=\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}$, with probability $1-\delta$, the approximate
solution $\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}$ is an $(\nu,\nu)$-approximate PCP and
\[
{\cal E}({\mathbf{x}}_{\mat R,k})\leq{\cal E}({\mathbf{x}}_{k})+\frac{(2\nu+\nu^{2})\TNormS{{\mathbf{f}}}}{n}
\]

Thus if, for example, $\ensuremath{{\bm{\mathrm{G}}}}$ is a \noun{CountSketch} matrix, then
the conditions are met when the number of rows in $\ensuremath{{\bm{\mathrm{G}}}}$ is
\[
\Omega\left(\frac{\sr{\ensuremath{{\bm{\mathrm{A}}}}}^{2}}{\gap{\ensuremath{{\bm{\mathrm{A}}}}}k^{2}\nu^{2}\delta}\right)
\]
 rows. In another example, if $\ensuremath{{\bm{\mathrm{G}}}}$ is a subgaussian map, then
the conditions are met when the number of rows in $\ensuremath{{\bm{\mathrm{G}}}}$ is
\[
\Omega\left(\frac{\sr{\ensuremath{{\bm{\mathrm{A}}}}}+\log(1/\delta)}{\gap{\ensuremath{{\bm{\mathrm{A}}}}}k^{2}\nu^{2}}\right)\,.
\]
\end{thm}

\begin{proof}
Due to Theorem~\ref{thm:structural} it suffices to show that that
$d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\leq\nu$.
Under the conditions of the Theorem, with probability of at least
$1-\delta$ we have $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}-\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}}\leq\epsilon\TNormS{\ensuremath{{\bm{\mathrm{A}}}}}$.
If that is indeed the case, $\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}$ has
rank at least $k$ since $\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}$ and $\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}$
are symmetric matrices and we know that $\sigma_{i}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\right)=\sigma_{i}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}\right)\pm\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}-\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}}$
(Weyl's Theorem and the fact that $\sigma_{k}^{2}>\epsilon\sigma_{1}^{2}$).
Furthermore, since $\nu>0$ we have $\epsilon<\gap{\ensuremath{{\bm{\mathrm{A}}}}}k$, and
Theorem~\ref{thm:sin-theta-1} implies
\begin{eqnarray*}
d_{2}(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}) & \leq & \frac{\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}-\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}}}{(\sigma_{k}^{2}-\sigma_{k+1}^{2})-\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}-\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}}}\\
 & \leq & \frac{\epsilon}{\gap{\ensuremath{{\bm{\mathrm{A}}}}}k-\epsilon}\\
 & \leq & \nu\,.
\end{eqnarray*}
Thus, we have shown that with probability $1-\delta$ we have $d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\leq\nu$,
as required.
\end{proof}
\begin{thm}
[Two Sided Sketching]\label{thm:left_right_sketching}Let $\nu,\delta\in(0,1/2)$
and denote
\[
\epsilon_{2}=\frac{\nu}{2(1+\nu/2)}\cdot\gap{\ensuremath{{\bm{\mathrm{A}}}}}k.
\]
Suppose that $\ensuremath{{\bm{\mathrm{G}}}}$ is sampled from a distribution that provides
a $(\epsilon_{2},\delta/2)$-approximate Gram matrix for $\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}$.
Denote
\[
\epsilon_{1}=\frac{\nu(1+\nu^{2}/4)^{-1/2}/2}{1+\nu(1+\nu^{2}/4)^{-1/2}/2}\cdot\gap{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}k
\]
and uppose that $\mat S$ is sampled from a distribution that provides
a $(\epsilon_{1},\delta/2)$-approximate Gram matrix for $\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}$.
Then for $\mat R=\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k}$ with probability
$1-\delta$ the approximate solution $\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}$ is an $(\nu,\nu)$-approximate
PCP and
\[
{\cal E}({\mathbf{x}}_{\mat R,k})\leq{\cal E}({\mathbf{x}}_{k})+\frac{(2\nu+\nu^{2})\TNormS{{\mathbf{f}}}}{n}
\]

Thus if, for example, $\mat S$ is a \noun{CountSketch} matrix and
$\ensuremath{{\bm{\mathrm{G}}}}$ is a subgaussian map, then the conditions hold when the numbers
of rows of $\mat S$ is
\[
\Omega\left(\frac{\sr{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}^{2}}{\gap{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}k^{2}\nu^{2}\delta}\right)
\]
 and the number of rows in $\ensuremath{{\bm{\mathrm{G}}}}$ is
\[
\Omega\left(\frac{\sr{\ensuremath{{\bm{\mathrm{A}}}}}+\log(1/\delta)}{\gap{\ensuremath{{\bm{\mathrm{A}}}}}k^{2}\nu^{2}}\right)\,.
\]
\end{thm}

\begin{proof}
Due to Theorem~\ref{thm:structural} it suffices to show that that
$d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\leq\nu$.

Under the conditions of the Theorem, with probability of at least
$1-\delta/2$ we have $\TNorm{(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}-(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}\leq\epsilon_{1}\TNormS{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}$,
and with probability of at least $1-\delta/2$ we have $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}-\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}}\leq\epsilon_{2}\TNormS{\ensuremath{{\bm{\mathrm{A}}}}}$.
Thus, both inequalities hold with probability of at least $1-\delta$.
If that is indeed the case, $\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}$ has
rank at least $k$ since $\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}$ and $\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}$
are symmetric matrices and we know that $\sigma_{i}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\right)=\sigma_{i}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}\right)\pm\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}-\ensuremath{{\bm{\mathrm{A}}}}\matA^{\textsc{T}}}$
(Weyl's Theorem and the fact that $\sigma_{k}^{2}>\epsilon_{2}\sigma_{1}^{2}$).
Moreover, $(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}$
has rank at least $k$ since $(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}$
and $\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}$ are symmetric matrices and we
know that $\sigma_{i}^{2}\left((\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\right)=\sigma_{i}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}\right)\pm\TNorm{(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}-\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{G}}}}\ensuremath{{\bm{\mathrm{A}}}}^{\textsc{T}}}$
(Weyl's Theorem and the fact that $\sigma_{k}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\right)>\epsilon_{1}\sigma_{1}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\right)$).
Since $\nu>0$ we have $\epsilon_{1}<\gap{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}k$ and
Theorem~\ref{thm:sin-theta-1} implies

\begin{eqnarray*}
d_{2}(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k}},\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k}) & \leq & \frac{\TNorm{(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}-(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}}{\sigma_{k}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\right)-\sigma_{k+1}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\right)-\TNorm{(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\mat S^{\textsc{T}}\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}-(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}}\\
 & \leq & \frac{\epsilon_{1}\TNormS{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}}{\sigma_{k}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\right)-\sigma_{k+1}^{2}\left(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\right)-\epsilon_{1}\TNormS{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}}\\
 & = & \frac{\epsilon_{1}}{\gap{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}k-\epsilon_{1}}\\
 & \leq & \nu(1+\nu^{2}/4)^{-1/2}/2
\end{eqnarray*}
From Lemma~\ref{lem:14} (with $\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}$) we get that $d_{2}(\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k}},\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k})\leq\frac{\sigma_{k+1}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})}{\sigma_{k}(\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})}(\nu/2)\leq\nu/2$.

We now bound
\begin{eqnarray*}
d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right) & \leq & d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k,k}\right)+d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\\
 & \leq & \nu
\end{eqnarray*}
where we similarly use Theorem \ref{thm:right_sketching} to bound
$d_{2}(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k})\leq\nu/2.$

Thus, we have shown that with probability $1-\delta$ we have $d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\leq\nu$,
as required.
\end{proof}

\subsection{Fast Approximate PCR/PCP}

A prototypical algorithm for approximate PCR/PCP is to compute ${\mathbf{x}}_{\mat R,k}$
with some choice of sketching-based $\mat R$. There are quite a few
design choices that need to be made in order to turn this prototypical
algorithm into a concrete algorithm, e.g. whether to use left, right
or two sided sketching to form $\mat R$, and which sketch transform
to use. There are various tradeoffs, e.g. using \noun{CountSketch}
results in faster sketching, but usually requires larger sketch sizes.
Furthermore, in computing ${\mathbf{x}}_{\mat R,k}$ there are also algorithmic
choices to be made with respect to choosing the order of matrix multiplications:
in computing $(\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{\textsc{+}}{\mathbf{b}}$ should
we first compute $\ensuremath{{\bm{\mathrm{A}}}}\mat R$ and then multiply by $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}$,
or vice versa? Likely, there is no one size fit all algorithm, and
different profiles of the input matrix (in particular, the size and
sparsity level) call for a different variant of the prototypical algorithm.

\begin{table}
\begin{centering}
{\footnotesize{}}%
\begin{tabular}{|c|c|c|c|}
\hline
 &  & {\footnotesize{}${\mathbf{x}}_{\mat R,k}$} & {\footnotesize{}CLS}\tabularnewline
\hline
\hline
{\footnotesize{}Left sketching} & {\footnotesize{}subgaussian $\mat S$} & {\footnotesize{}$s_{1}\cdot\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+s_{1}d\min(s_{1},d)+nk^{2}$} & {\footnotesize{}N/A}\tabularnewline
\cline{2-4}
{\footnotesize{}$\mat R=\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}},k}$} & \noun{\footnotesize{}CountSketch}{\footnotesize{} $\mat S$} & {\footnotesize{}$k\cdot\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+s_{2}d\min(s_{2},d)+nk^{2}$} & {\footnotesize{}N/A}\tabularnewline
\hline
{\footnotesize{}Right sketching} & {\footnotesize{}subgaussian $\ensuremath{{\bm{\mathrm{G}}}}$} & {\footnotesize{}$t_{1}\cdot\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+nt_{1}\min(n,t_{1})+t_{1}k\min(n,d)$} & {\footnotesize{}$t_{1}\cdot\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+nt_{1}\min(n,t_{1})$}\tabularnewline
\cline{2-4}
{\footnotesize{}$\mat R=\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}$} & \noun{\footnotesize{}CountSketch}{\footnotesize{} $\ensuremath{{\bm{\mathrm{G}}}}$} & {\footnotesize{}$\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+nt_{2}\min(n,t_{2})+t_{2}k\min(n,d)$} & {\footnotesize{}$\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+nt_{2}\min(n,t_{2})$}\tabularnewline
\hline
{\footnotesize{}Two sided } & \noun{\footnotesize{}CountSketch} & {\footnotesize{}$\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+s_{2}k^{2}+k\min(nt_{2},\nnz{\ensuremath{{\bm{\mathrm{A}}}}})+nk^{2}$} & {\footnotesize{}N/A}\tabularnewline
{\footnotesize{}$\mat R=\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k}$} & {\footnotesize{}$\ensuremath{{\bm{\mathrm{G}}}}$ and $\mat S$} &  & \tabularnewline
\hline
\end{tabular}
\par\end{centering}{\footnotesize \par}
\caption{\label{tab:complex}Computational complexity of computing $\protect{\mathbf{x}}_{\protect\mat R,k}$
and CLS for various options $\protect\mat R$. For brevity, we omit
the $O$$()$ from the notation. }
\end{table}

Table~\ref{tab:complex} summarizes the running time complexity of
several design options. In order to better make sense between these
different choices, we first summarize the running time complexity
of various design choices using the optimal implementation (from an
asymptotic running-time complexity perspective). To make the discussion
manageable, we consider only subgaussian maps and \noun{CountSketch}.
Furthermore, for the sake of the analysis, we make some assumptions
and adopt some notational conventions. First, we assume that computing
$\ensuremath{{\bm{\mathrm{B}}}}^{\textsc{+}}{\mathbf{c}}$ for some $\ensuremath{{\bm{\mathrm{B}}}}\in\mathbb{R}^{m\times n}$ and ${\mathbf{c}}$ is
done via straightforward methods based on QR or SVD factorizations,
and as such takes $O(mn\min(m,n))$. We consider using fast sketch-based
approximate least squares algorithms in the next subsection. Next,
we let the sketch sizes be parameters in the complexity. In the discussion,
we use our theoretical results to deduce reasonable assumptions on
how these parameters are set, and thus to reason about the final complexity
of sketched PCR/PCP. We denote the number of rows in the left sketch
matrix $\mat S$ by $s_{1}$ for a subgaussian map, and $s_{2}$ for
\noun{CountSketch}. We denote the number of rows in the left sketch
matrix $\ensuremath{{\bm{\mathrm{G}}}}$ by $t_{1}$ for a subgaussian map, and $t_{2}$ for
\noun{CountSketch}. Finally, we assume $\nnz{\ensuremath{{\bm{\mathrm{A}}}}}\geq\max(n,d)$,
and that all sketch sizes are greater than $k$.

Table~\ref{tab:complex} also lists, where relevant, the complexity
of the CLS solution ${\mathbf{x}}_{\mat R}$.

\paragraph{Discussion.}

We first compare the computational complexity of CLS to the computational
complexity of our proposed right sketching algorithm. For both choices
of $\ensuremath{{\bm{\mathrm{G}}}}$ we have for sketched PCP an additional term of $O(tk\min(n,d))$.
However, close inspection reveals that this term is dominated by the
term $O(nt\min(n,t))$. Thus our proposed algorithm has the same asymptotic
complexity as CLS for the same sketch size. However, our algorithm
does not mix regularization and compression and comes with stronger
theoretical guarantees.

Next, in order to compare subgaussian maps to \noun{CountSketch},
we first make some simplified assumptions on the required approximation
quality $\nu$, the relative eigengap $\gap{\ensuremath{{\bm{\mathrm{A}}}}}k$, and the rank
parameter $k$: $\nu$ is fixed, $\gap{\ensuremath{{\bm{\mathrm{A}}}}}k$ is bounded from below
by a constant, and we have $k=O(\sr{\ensuremath{{\bm{\mathrm{A}}}}})$. The first assumptions
is justified if we are satisfied with fixed sub-optimality in the
objective (optimization perspective), or a small constant multiplicative
increase in excess risk if left sketching is used, or $n$ is fixed
(statistical perspective). The first assumption is somewhat less justified
from a statistical point of view when $n\to\infty$ and right sketching
is used. The rationale behind the second assumption is that the PCR/PCP
problem is in a sense ill-posed if there is a tiny eigengap. The third
assumption is motivated by the fact that the stable rank is a measure
of the number of large singular values, which are typically singular
values that correspond to the signal rather than noise. With these
assumptions, our theoretical results establish that $s_{1},t_{1}=O(k)$
and $s_{2},t_{2}=O(k^{2})$ suffice. It is important to stress that
we make these assumptions only for the sake of comparing the different
sketching options, and we do not claim that these assumptions always
hold, or that our proposed algorithms work only when these assumptions
hold.

For left sketching, with these assumptions, we have complexity of
$O(k\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+k^{2}\max(n,d))$ for subgaussian maps and $O(k\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+dk^{2}\min(k^{2},d)+nk^{2})$
for \noun{CountSketch}. Clearly, better asymptotic complexity is achieved
with subgaussian maps. For right sketching, with these assumptions,
we have complexity of $O(k\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+nk^{2})$ for subgaussian sketch
and $O(\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+nk^{2}\min(n,k^{2}))$ for \noun{CountSketch}.
The complexity in terms of the input sparsity $\nnz{\ensuremath{{\bm{\mathrm{A}}}}}$\noun{,}
which is arguably the dominant term, is better for \noun{CountSketch}.
For two sided sketching, we have complexity $O(\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+k^{4}+k\min(nk^{2},\nnz{\ensuremath{{\bm{\mathrm{A}}}}})+nk^{2}$).

If $n\gg d$ and $\nnz{\ensuremath{{\bm{\mathrm{A}}}}}=O(n)$ (sparse input matrix, and constant
amount of non zero features per data point), left sketching gives
better asymptotic complexity. If $n\gg d$ and $\nnz{\ensuremath{{\bm{\mathrm{A}}}}}=nd$ (full
data matrix), left sketching has better complexity unless $d\gg k^{3}$.
Furthermore, left sketching gives stronger theoretical guarantees.\emph{
Thus, for $n\gg d$} \emph{we advocate the use of left sketching}.
If $d\gg n$ and $\nnz{\ensuremath{{\bm{\mathrm{A}}}}}=O(d)$ (sparse input matrix), right
sketching with a subgaussian maps always has better complexity than
left sketching, and potentially (but not always) right sketching with
\noun{CountSketch} has even better complexity. If $d\gg n$ and $\nnz{\ensuremath{{\bm{\mathrm{A}}}}}=nd$
(full data matrix), right sketching with subgaussian maps has the
same complexity as left sketching, and potentially (but not always)
right sketching with \noun{CountSketch} has even better complexity
(if $d$ is sufficiently larger than $n$). \emph{Thus, for $d\gg n$
we advocate the use of right sketching}. If $n\approx d$ (both very
large), and $\nnz{\ensuremath{{\bm{\mathrm{A}}}}}=n$ then it is possible to have $O(nk^{2})$
with all three options (left, right and two sided), as long as $k^{2}\leq n$.
A similar conclusion is achieved if $n\approx d$ and $\nnz{\ensuremath{{\bm{\mathrm{A}}}}}=nd$,
but if $k^{2}\ll n$ then two sided sketch is better.

\subsection{\label{subsec:input-sparsity}Input Sparsity Approximate PCP}

In this section, we propose an input-sparsity algorithm for approximate
PCP. In 'input-sparsity algorithm', we mean an algorithm whose running
time is $O(\nnz{\ensuremath{{\bm{\mathrm{A}}}}}\log(d/\epsilon)+\poly{k,s,t,\log(1/\epsilon})$,
where $\epsilon$ is some accuracy parameter (see formal theorem statement).

\begin{algorithm}
\begin{algorithmic}[1]

\STATE \textbf{Input: $\ensuremath{{\bm{\mathrm{A}}}}\in\mathbb{R}^{n\times d}$},\textbf{ ${\mathbf{b}}\in\mathbb{R}^{n}$},\textbf{
}$k\leq\min(n,d)$, $s,t\geq k$, $\epsilon\in(0,1)$

\STATE

\STATE Generate two \noun{CountSketch} matrices $\mat S\in\mathbb{R}^{s\times n}$
and $\ensuremath{{\bm{\mathrm{G}}}}\in\mathbb{R}^{t\times n}$.

\STATE Remove from $\ensuremath{{\bm{\mathrm{G}}}}$ any row that is zero.

\STATE $\ensuremath{{\bm{\mathrm{C}}}}\gets\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}$.

\STATE $\ensuremath{{\bm{\mathrm{D}}}}\gets\mat S\ensuremath{{\bm{\mathrm{C}}}}$.

\STATE Compute $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{D}}}},k}$, the $k$ dominant right invariant
space of $\ensuremath{{\bm{\mathrm{G}}}}$ (via SVD).

\STATE For the analysis (no need to compute): $\mat R=\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{D}}}},k}$.

\STATE Solve $\min_{\gamma}\TNorm{\ensuremath{{\bm{\mathrm{C}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{D}}}},k}\gamma-{\mathbf{b}}}$
to $\epsilon/d$ accuracy using input-sparsity least squares regression
(see \cite[section 7.7]{CW17}). (Do not compute $\ensuremath{{\bm{\mathrm{C}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{D}}}},k}$.
In each iteration, multiplying a vector by $\ensuremath{{\bm{\mathrm{C}}}}\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{D}}}},k}$
is performed by first multiplying by $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{D}}}},k}$ and then by
$\ensuremath{{\bm{\mathrm{C}}}}$.)

\STATE Return ${\mathbf{y}}\gets\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}(\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{D}}}},k}\tilde{\gamma})$,
where $\tilde{\gamma}$ is the output of the previous step.

\end{algorithmic}

\caption{\label{alg:input_sparsity}Input Sparsity Approximate PCP}
\end{algorithm}

The basic idea is to use two sided sketching, with an additional modification
of using input sparsity algorithms to approximate $(\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{\textsc{+}}{\mathbf{b}}=\arg\min_{\gamma}\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\mat R\gamma-{\mathbf{b}}}$.
Specifically, we propose to use the algorithm recently suggested by
Clarkson and Woodruff~\cite{CW17}. A pseudo-code description of
our input sparsity approximate PCP algorithm is listed in Algorithm~\ref{alg:input_sparsity}.
We have the following statement about the algorithm.
\begin{thm}
Run Algorithm \ref{alg:input_sparsity} with $\epsilon,s,t,k$ as
parameters. Under exact arithmetic\footnote{The results are likely too optimistic for inexact arithmetic. We leave
the numerical analysis to future work. }, after
\[
O(\nnz{\ensuremath{{\bm{\mathrm{A}}}}}\log(d/\epsilon)+\log(d/\epsilon)tk+sk^{2}+t^{3}k+k^{3}\log^{2}k)
\]
operations, with probability $2/3$, the algorithm will return a ${\mathbf{y}}$
such that
\[
\TNormS{{\mathbf{y}}-{\mathbf{x}}_{\mat R}}\leq\epsilon\TNormS{{\mathbf{x}}_{\mat R}}\,.
\]
\end{thm}

\begin{proof}
Denote $\ensuremath{{\bm{\mathrm{B}}}}=\ensuremath{{\bm{\mathrm{A}}}}\mat R$, and consider using the iterative method
described in~\cite[section 7.7]{CW17} to approximately solve $\min_{\gamma}\TNorm{\ensuremath{{\bm{\mathrm{B}}}}\gamma-{\mathbf{b}}}$.
Denote the optimal solution by $\gamma_{\mat R}$, and the solution
that our algorithm found by $\tilde{\gamma}$. Theorem 7.14 in \cite{CW17}
states that after the $O(\log(d/\epsilon))$ iterations the algorithm
would have returned $\tilde{\gamma}$ such that
\begin{equation}
\TNormS{\ensuremath{{\bm{\mathrm{B}}}}\mat Z(\tilde{\gamma}-\gamma_{\mat R})}\leq(\epsilon/d)\TNormS{\ensuremath{{\bm{\mathrm{B}}}}\mat Z\gamma_{\mat R}}\label{eq:gamma_err}
\end{equation}
for some invertible $\mat Z$ found by the algorithm. Furthermore,
$\kappa(\ensuremath{{\bm{\mathrm{B}}}}\mat Z)=O(1)$ where $\kappa(\cdot)$ is the condition
number (ratio between the largest singular value and smallest). Eq.~(\ref{eq:gamma_err})
implies that $\TNormS{\tilde{\gamma}-\gamma_{\mat R}}\leq\kappa(\ensuremath{{\bm{\mathrm{B}}}}\mat Z)^{2}(\epsilon/d)\TNormS{\gamma_{\mat R}}=O(\epsilon/d)\TNormS{\gamma_{\mat R}}$.
Now, noticing that ${\mathbf{x}}_{\mat R}=\mat R\gamma_{\mat R}$ and ${\mathbf{y}}=\mat R\tilde{\gamma}$,
we find that
\[
\TNormS{{\mathbf{y}}-{\mathbf{x}}_{\mat R}}\leq O(\epsilon/d)\kappa(\mat R)^{2}\TNormS{{\mathbf{x}}_{\mat R}}\,.
\]
We now need to bound $\kappa(\mat R)=\kappa(\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}}\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}},k})=\kappa(\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})$
where $\ensuremath{{\bm{\mathrm{G}}}}$ is a \noun{CountSketch} matrix. Since $\ensuremath{{\bm{\mathrm{G}}}}$ has
a single non zero in each column, then $\TNormS{\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}\leq\FNormS{\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}}\leq d$.
Furthermore, since we removed zero column from $\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}$, for
any ${\mathbf{x}}$ the vector $\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}}{\mathbf{x}}$ has in one of its coordinates
any coordinate of ${\mathbf{x}}$, so $\sigma_{\min}(\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})\geq1$. So
we found that $\kappa(\ensuremath{{\bm{\mathrm{G}}}}^{\textsc{T}})^{2}\leq d$. We conclude that
\[
\TNormS{{\mathbf{y}}-{\mathbf{x}}_{\mat R}}\leq O(\epsilon)\TNormS{{\mathbf{x}}_{\mat R}}\,.
\]
Adjusting $\epsilon$ to compensate for the constants completes the
proof.
\end{proof}

\section{\label{sec:extensions}Extensions}

\subsection{Streaming Algorithm}

We now consider computing an approximate PCR/PCP in the streaming
model. We consider a one-pass row-insertion streaming model, in which
the rows of $\ensuremath{{\bm{\mathrm{A}}}}$, ${\mathbf{a}}_{1},{\mathbf{a}}_{2},\dots,{\mathbf{a}}_{n}$, and the corresponding
entries in ${\mathbf{b}}$, $b_{1},b_{2},\dots,b_{n}$, are presented one-by-one
and once only (i.e., in a stream). The goal is to use $o(n)$ memory
(so $\ensuremath{{\bm{\mathrm{A}}}}$ cannot be stored in memory). The relevant resources to
be bounded for numerical linear algebra in the streaming model are
storage, update time (time spent per row), and final computation time
(at the end of the stream)~\cite{CW09}. Our goal is to bound these
by $O(\poly d)$ .

Our proposed streaming algorithm for approximate PCP uses left sketching.
It is easy to verify that if $\mat S$ is a subgaussian map or \noun{CountSketch,
}then $\mat R=\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}},k}$ can be computed in the streaming
model: one has to update $\mat S\ensuremath{{\bm{\mathrm{A}}}}$ as new rows are presented ($O(d)$
update for \noun{CountSketch}, and $O(sd)$ for subgaussian map),
and once the final row has been presented, factorizing $\mat S\ensuremath{{\bm{\mathrm{A}}}}$
and extracting $\mat R$ can be done in $O(sd\min(s,d))$ which is
polynomial in $d$ if $s$ is polynomial in $d$. However, to compute
${\mathbf{x}}_{\mat R}$ one has to compute $(\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{\textsc{+}}{\mathbf{b}}$, and storing
$\ensuremath{{\bm{\mathrm{A}}}}\mat R$ in memory requires $\Omega(n)$ memory. To circumvent
this issue we propose to introduce another sketching matrix $\mat T$,
and approximate $(\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{\textsc{+}}{\mathbf{b}}$ via $(\mat T\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{\textsc{+}}{\mathbf{b}}$
. Thus, for $\mat R=\ensuremath{{\bm{\mathrm{V}}}}_{\mat S\ensuremath{{\bm{\mathrm{A}}}},k}$ we approximate ${\mathbf{x}}_{\mat R}$
by $\tilde{{\mathbf{x}}}_{\mat R}=\mat R(\mat T\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{\textsc{+}}{\mathbf{b}}$. It is easy
to verify that $\tilde{{\mathbf{x}}}_{\mat R}$ can be computed in the streaming
model (by forming and updating $\mat T\ensuremath{{\bm{\mathrm{A}}}}$ while computing $\mat R$).

More generally, for \emph{any} $\mat R$ which can be computed in the
streaming model, we can also compute in the streaming model the following
approximation of ${\mathbf{x}}_{\mat R,k}$:

\[
\tilde{{\mathbf{x}}}_{\mat R,k}:=\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}(\mat T\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}\mat T{\mathbf{b}}\,.
\]
The next theorem, establishes conditions on $\mat T$ that guarantee
that $\tilde{{\mathbf{x}}}_{\mat R,k}$ is a an approximate PCR/PCP.
\begin{thm}
\label{thm:sketch-stream-structure-impl}Suppose that $\mat R\in\mathbb{R}^{d\times s}$
with $s\ge k$. Assume that $\nu\in(0,1)$. Suppose that $\mat T$
provides a $O(\nu)$-distortion subspace embedding for $\range{\left[\begin{array}{ccc}
\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k} & \ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k} & {\mathbf{b}}\end{array}\right]}$ that is
\[
\TNormS{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}{\mathbf{x}}_{1}+\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}{\mathbf{x}}_{2}+{\mathbf{b}} x_{3}}=(1\pm O(\nu))(\TNormS{{\mathbf{x}}_{1}}+\TNormS{{\mathbf{x}}_{2}}+x_{3}^{3})
\]
for every ${\mathbf{x}}_{1},{\mathbf{x}}_{2}\in\mathbb{R}^{k}$ and $x_{3}\in\mathbb{R}$. Then,
\begin{enumerate}
\item If $d_{2}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\mat A\mat R,k},\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\right)\leq\nu$
then $\ensuremath{{\bm{\mathrm{A}}}}\tilde{{\mathbf{x}}}_{\mat R,k}$ is an $(O(\nu),O(\nu))$-approximate
PCP.
\item If $s=k$ and $\mat R$ has orthonormal columns (i.e., $\mat R^{\textsc{T}}\mat R=\mat I_{k})$
and $d_{2}\left(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k}\right)\leq\nu(1+\nu^{2})^{-1/2}$
then $\tilde{{\mathbf{x}}}_{\mat R}$ is an $(O(\nu),\,O(\nu/\sigma_{k}))$-approximate
PCR.
\end{enumerate}
The subspace embedding conditions on $\mat T$ are met with probability
of at least $1-\delta$ if, for example, $\mat T$ is a \noun{CountSketch}
matrix with $O(k^{2}/\nu^{2}\delta)$ rows.
\end{thm}

\begin{proof}
We need to show both the additive error bounds on the objective function,
and the error bound on the constraints. We start with the additive
error bounds on the objective function for both for PCP (first part
of the theorem) and for PCR (second part of the theorem). The lower
bound on $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\tilde{{\mathbf{x}}}_{\mat R,k}-{\mathbf{b}}}$ follows immediately
from the fact that $\tilde{{\mathbf{x}}}_{\mat R,k}\in\range{\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}$
and the fact that ${\mathbf{x}}_{\mat R,k}$ is a minimizer of $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}-{\mathbf{b}}}$
subject to ${\mathbf{x}}\in\range{\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}}$. For the upper
bound, we observe
\begin{eqnarray*}
\TNorm{\ensuremath{{\bm{\mathrm{A}}}}\tilde{{\mathbf{x}}}_{\mat R,k}-{\mathbf{b}}} & \leq & \left(1+O(\nu)\right)\TNorm{\mat T\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}(\mat T\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}\mat T{\mathbf{b}}-\mat T{\mathbf{b}}}\\
 & \leq & \left(1+O(\nu)\right)\TNorm{\mat T\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}(\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}{\mathbf{b}}-\mat T{\mathbf{b}}}\\
 & \leq & \left(1+O(\nu)\right)\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}-{\mathbf{b}}}\\
 & \leq & \TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}-{\mathbf{b}}}+O(\nu)\TNorm{{\mathbf{b}}}
\end{eqnarray*}
where in the first and third inequality we used the fact that $\mat T$
provides a subspace embedding for $\range{[\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\,{\mathbf{b}}]}$
and in the second inequality we used the fact that $(\mat T\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}\mat T{\mathbf{b}}$
is a minimizer of $\TNorm{\mat T\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}{\mathbf{x}}-\mat T{\mathbf{b}}}$.
Bounds on $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{\mat R,k}-{\mathbf{b}}}$ (Theorem~\ref{thm:structural})
now imply the additive bound.

We now bound the constraint for the PCR guarantee (second part of
the theorem). Let
\[
\ensuremath{{\bm{\mathrm{C}}}}=(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}((\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{\textsc{T}})^{+}\,.
\]
 Since $(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{\textsc{T}}$ and $(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}$
have the same row space, and $\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}$ has more
rows than columns, $\ensuremath{{\bm{\mathrm{C}}}}$ is non-singular and we have $\ensuremath{{\bm{\mathrm{C}}}}(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{\textsc{T}}=(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}$.
Since $\mat T$ provides a subspace embedding for $\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}$,
all the singular values of $\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}$ belong to
the interval $[1-O(\nu),1+O(\nu)]$. We conclude that $\TNorm{\ensuremath{{\bm{\mathrm{C}}}}-\mat I_{k}}\leq O(\nu)$.
We also have $(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}=\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}^{-1}(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}$
since $\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}$ has linearly independent columns
(since it provides a subspace embedding), and $\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}$
has all linearly independent rows. Thus,

\begin{eqnarray*}
\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}\tilde{{\mathbf{x}}}_{\mat R,k}} & = & \TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{A}}}}\tilde{{\mathbf{x}}}_{\mat R,k}-\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k+}^{\textsc{T}}\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}^{\textsc{T}}\mat T^{\textsc{T}}\mat T{\mathbf{b}}}\\
 & \leq & \TNorm{\ensuremath{{\bm{\mathrm{A}}}}\tilde{{\mathbf{x}}}_{\mat R,k}-\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}^{\textsc{T}}\mat T^{\textsc{T}}\mat T{\mathbf{b}}}\\
 & = & \TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}\mat T{\mathbf{b}}-\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}^{\textsc{T}}\mat T^{\textsc{T}}\mat T{\mathbf{b}}}\\
 & = & \TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\ensuremath{{\bm{\mathrm{C}}}}(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{\textsc{T}}\mat T{\mathbf{b}}-\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}^{\textsc{T}}\mat T^{\textsc{T}}\mat T{\mathbf{b}}}\\
 & \leq & \left(1+O(\nu)\right)\cdot\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\ensuremath{{\bm{\mathrm{C}}}}\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}^{\textsc{T}}\mat T^{\textsc{T}}-\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}^{\textsc{T}}\mat T^{\textsc{T}}}\cdot\TNorm{{\mathbf{b}}}\\
 & \leq & \left(1+O(\nu)\right)^{2}\cdot\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\ensuremath{{\bm{\mathrm{C}}}}\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}^{\textsc{T}}-\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}^{\textsc{T}}}\TNorm{{\mathbf{b}}}\\
 & \leq & \left(1+O(\nu)\right)\cdot\left(\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\left(\ensuremath{{\bm{\mathrm{C}}}}-\mat I_{k}\right)\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}^{\textsc{T}}}+\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}^{\textsc{T}}-\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}\ensuremath{{\bm{\mathrm{U}}}}_{\mat A,k}^{\textsc{T}}}\right)\cdot\TNorm{{\mathbf{b}}}\\
 & \leq & \left(1+O(\nu)\right)\cdot\left(\TNorm{\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\left(\ensuremath{{\bm{\mathrm{C}}}}-\mat I_{k}\right)\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}^{\textsc{T}}}+\nu\right)\cdot\TNorm{{\mathbf{b}}}\\
 & = & \left(1+O(\nu)\right)\cdot\left(\TNorm{\left(\ensuremath{{\bm{\mathrm{C}}}}-\mat I_{k}\right)}+\nu\right)\cdot\TNorm{{\mathbf{b}}}\\
 & \leq & \left(1+O(\nu)\right)\cdot\left(O(\nu)+\nu\right)\cdot\TNorm{{\mathbf{b}}}\\
 & = & O(\nu)\cdot\TNorm{{\mathbf{b}}}
\end{eqnarray*}

We now bound the constraint for the PCR guarantee (second part of
the theorem). To that end, and observe:

\begin{eqnarray*}
\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\tilde{{\mathbf{x}}}_{\mat R,k}} & \leq & \TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}(\mat T\ensuremath{{\bm{\mathrm{A}}}}\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}\mat T{\mathbf{b}}}\\
 & \leq & \TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R}\cdot\TNorm{(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k})^{+}\mat T{\mathbf{b}}}\\
 & \leq & \frac{\nu\cdot\left(1+O(\nu)\right)\cdot\TNorm{{\mathbf{b}}}}{\sigma_{\min}\left(\mat T\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\right)}\\
 & \leq & \frac{\nu\cdot\left(1+O(\nu)\right)\TNorm{{\mathbf{b}}}}{\left(1-O(\nu\right))\sigma_{\min}\left(\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\mat{\Sigma}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k}\right)}\\
 & \leq & \frac{O(\nu)\cdot\TNorm{{\mathbf{b}}}}{\sigma_{\min}\left(\ensuremath{{\bm{\mathrm{A}}}}\mat R\right)}\\
 & \leq & \frac{O(\nu)}{\sigma_{k}}\cdot\TNorm{{\mathbf{b}}}
\end{eqnarray*}
where we used the fact that $\mat T$ provides a subspace embedding
for $\range{\left[\begin{array}{cc}
\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{A}}}}\mat R,k} & {\mathbf{b}}\end{array}\right]}$, and used Lemma~\ref{lem:11} to bound $\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\mat A,k+}^{\textsc{T}}\mat R}$
and $\TNorm{(\ensuremath{{\bm{\mathrm{A}}}}\mat R)^{+}}$.
\end{proof}

\subsection{Approximate Kernel PCR}

For simplicity, we consider only the homogeneous polynomial kernel
${\cal K}({\mathbf{x}},{\mathbf{z}})=({\mathbf{x}}^{\textsc{T}}{\mathbf{z}})^{q}$. The results trivially extend to
the non-homogeneous polynomial kernel ${\cal K}_{n}({\mathbf{x}},{\mathbf{z}})=({\mathbf{x}}^{\textsc{T}}{\mathbf{z}}+c)^{q}$
by adding a single feature to each data point. We leave to future
work the development of similar techniques for other kernels (e.g.
Gaussian kernel).

Let $\phi:\mathbb{R}^{d}\to\mathbb{R}^{d^{q}}$ be the function that maps a vector
${\mathbf{z}}=(z_{1},\dots,z_{d})$ to the set of monomials formed by multiplying
$q$ entries of ${\mathbf{z}}$, i.e. $\phi({\mathbf{z}})=(z_{i_{1}}z_{i_{2}}\cdots z_{i_{q}})_{i_{1},\dots,i_{q}\in\{1,\dots,d\}}$.
For a data matrix $\ensuremath{{\bm{\mathrm{A}}}}\in\mathbb{R}^{d}$ and a response vector ${\mathbf{b}}\in\mathbb{R}^{n}$,
let $\Phi\in\mathbb{R}^{n\times d^{q}}$ be the matrix obtained by applying
$\phi$ to the rows of $\ensuremath{{\bm{\mathrm{A}}}}$, and consider computing the rank $k$
PCR solution $\Phi$ and ${\mathbf{b}}$, which we denote by ${\mathbf{x}}_{{\cal K},k}$.
The corresponding prediction function is $f_{{\cal K},k}({\mathbf{z}})\coloneqq\phi({\mathbf{z}})^{\textsc{T}}{\mathbf{x}}_{{\cal K},k}$.
While ${\mathbf{x}}_{{\cal K},k}$ is likely a huge vector (since ${\mathbf{x}}_{{\cal K},k}\in\mathbb{R}^{d^{q}}$),
and thus expensive to compute, in kernel PCR we are primarily interested
in having an efficient method to compute $f_{{\cal K},k}({\mathbf{z}})$ given
a 'new' ${\mathbf{z}}$. We can accomplish this via the kernel trick, as we
now show.

We assume that $\Phi$ has full row rank (this holds if all data points
are different). Let ${\mathbf{a}}_{1},\dots,{\mathbf{a}}_{n}$ be the rows of $\ensuremath{{\bm{\mathrm{A}}}}$.
As usual with PCR, we have ${\mathbf{x}}_{{\cal K},k}=\ensuremath{{\bm{\mathrm{V}}}}_{\Phi,k}\Sigma_{\Phi,k}^{-1}\ensuremath{{\bm{\mathrm{U}}}}_{\Phi,k}^{\textsc{T}}{\mathbf{b}}$.
Since $\ensuremath{{\bm{\mathrm{V}}}}_{\Phi,k}=\Phi^{\textsc{T}}\ensuremath{{\bm{\mathrm{U}}}}_{\Phi,k}\Sigma_{\Phi,k}^{-1}$
we have
\begin{equation}
f_{{\cal K},k}({\mathbf{z}})=\phi({\mathbf{z}})^{\textsc{T}}\Phi^{\textsc{T}}\ensuremath{{\bm{\mathrm{U}}}}_{\Phi,k}\Sigma_{\Phi,k}^{-2}\ensuremath{{\bm{\mathrm{U}}}}_{\Phi,k}{\mathbf{b}}=({\cal K}({\mathbf{z}},{\mathbf{a}}_{1})\,\cdots\,{\cal K}({\mathbf{z}},{\mathbf{a}}_{n}))\alpha_{{\cal K},k}\label{eq:kpcr}
\end{equation}
where $\alpha_{{\cal K},k}\coloneqq\ensuremath{{\bm{\mathrm{U}}}}_{\Phi,k}\Sigma_{\Phi,k}^{-2}\ensuremath{{\bm{\mathrm{U}}}}_{\Phi,k}^{\textsc{T}}{\mathbf{b}}$.
In the above, we used the fact that for any ${\mathbf{x}}$ and ${\mathbf{z}}$ we have
$\phi({\mathbf{x}})^{\textsc{T}}\phi({\mathbf{z}})=({\mathbf{x}}^{\textsc{T}}{\mathbf{z}})^{q}={\cal K}({\mathbf{x}},{\mathbf{z}})$. Let $\ensuremath{{\bm{\mathrm{K}}}}\in\mathbb{R}^{n\times n}$
be the\emph{ kernel matrix} (also called \emph{Gram matrix}) defined
by $\ensuremath{{\bm{\mathrm{K}}}}_{ij}={\cal K}({\mathbf{a}}_{i},{\mathbf{a}}_{j}$). It is easy to verify that
$\ensuremath{{\bm{\mathrm{K}}}}=\Phi\Phi^{\textsc{T}}$, so we can compute $\ensuremath{{\bm{\mathrm{K}}}}$ in $O(n^{2}(d+q))$
(and without forming $\Phi$, which is a huge matrix). We also have
$\ensuremath{{\bm{\mathrm{K}}}}=\ensuremath{{\bm{\mathrm{U}}}}_{\Phi}\Sigma_{\Phi}^{2}\ensuremath{{\bm{\mathrm{U}}}}_{\Phi}^{\textsc{T}}$ so $\alpha_{k}=\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{K}}}},k}\Sigma_{K,k}^{-1}\ensuremath{{\bm{\mathrm{U}}}}_{\ensuremath{{\bm{\mathrm{K}}}},k}^{\textsc{T}}{\mathbf{b}}$.
Thus, we can compute $\alpha_{{\cal K},k}$ in $O(n^{2}(d+q+n))$
time. Once we have computed $\alpha_{k}$, using~(\ref{eq:kpcr})
we can compute $f_{{\cal K},k}({\mathbf{z}})$ for any ${\mathbf{z}}$ in $O(ndq)$ time.

In order to compute an approximate kernel PCR, we introduce a right
sketching matrix $\mat R\in\mathbb{R}^{d^{q}\times t}$. Such a matrix $\mat R$
is frequently referred to, in the context of kernel learning, as a
randomized feature map. We use the \noun{TensorSketch} feature map~\cite{Pagh13,PP13}.
The feature map is defined as follows. We first randomly generate
$q$ $3$-wise independent hash functions $h_{1},\dots,h_{q}\in\{1,\dots,d\}\to\{1,\dots,t\}$
and $q$ $4$-wise independent sign functions $g_{1},\dots,g_{q}:\{1,\dots,d\}\to\{-1,+1\}$.
Next, we define $H:\{1,\dots,d\}^{q}\to\{1,\dots,t\}$ and $G:\{1,\dots,t\}^{q}\to\{-1,+1\}$:
\[
H(i_{1},\dots,i_{q})\coloneqq h_{1}(i_{1})+\dots+h_{q}(i_{q})\,\mod\,t
\]
\[
G(i_{1},\dots,i_{q})=g_{1}(i_{1})\cdot g_{2}(i_{2})\cdot\dots\cdot g_{q}(i_{q})
\]
 To define $\mat R$, we index the rows of $\mat R$ by $\{1,\dots,d\}^{q}$
and set row $(i_{1},\dots,i_{q})$ to be equal to $G(i_{1},\dots,i_{q})\cdot{\mathbf{e}}_{H(i_{1},\dots,i_{q})}$,
where ${\mathbf{e}}_{j}$ denote the $j$th identity vector. A crucial observation
that makes \noun{TensorSketch} useful, is that via the representation
using $h_{1},\dots,h_{q}$ and $g_{1},\dots,g_{q}$ we can compute
$\mat R^{\textsc{T}}\phi({\mathbf{z}})$ in time $O(q(\nnz{{\mathbf{z}}}+t\log t))$ (see Pagh~\cite{Pagh13}
for details). Thus, we can compute $\Phi\mat R$ in time $O(q(\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+nt\log t))$.

Consider right sketching PCR on $\Phi$ and $k$ with a \noun{TensorSketch}
$\mat R$ as the sketching matrix. The approximate solution is
\[
{\mathbf{x}}_{{\cal K},\mat R,k}\coloneqq\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\Phi\mat R,k}(\Phi\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\Phi\mat R,k})^{\textsc{+}}{\mathbf{b}}=\mat R\gamma_{{\cal K},\mat R,k}
\]
 where $\gamma_{{\cal K},\mat R,k}\coloneqq\ensuremath{{\bm{\mathrm{V}}}}_{\Phi\mat R,k}(\Phi\mat R\ensuremath{{\bm{\mathrm{V}}}}_{\Phi\mat R,k})^{\textsc{+}}{\mathbf{b}}$.
We can compute $\gamma_{\mat R,k}$ in $O(q(\nnz{\ensuremath{{\bm{\mathrm{A}}}}}+nt\log t)+nt^{2})$
time. The predication function is
\[
f_{{\cal K},\mat R,k}({\mathbf{z}})\coloneqq\phi({\mathbf{z}})^{\textsc{T}}{\mathbf{x}}_{{\cal K},\mat R,k}=(\mat R^{\textsc{T}}\phi({\mathbf{z}}))^{\textsc{T}}\gamma_{{\cal K},\mat R,k}
\]
so once we have $\gamma_{{\cal K},\mat R,k}$ we can compute $f_{{\cal K},\mat R,k}({\mathbf{z}})$
in $O(q(\nnz{{\mathbf{z}}}+t\log t))$ time. Thus, the method is attractive
from a computational complexity point of view if $t\ll n$ or $d\gg n$
and $d\gg t$. The following theorem bound the excess risk of ${\mathbf{x}}_{{\cal K},\mat R,k}$.
\begin{thm}
Let $(\nu,\delta)\in(0,1/2)$. Let $\lambda_{1}\geq\dots\geq\lambda_{n}$
denote the eigenvalues of $\ensuremath{{\bm{\mathrm{K}}}}$. If $\mat R$ is a \noun{TensorSketch}
matrix with
\[
t=\Omega\left(\frac{3^{q}\Trace{\ensuremath{{\bm{\mathrm{K}}}}}^{2}}{(\lambda_{k}-\lambda_{k+1})^{2}\nu^{2}\delta}\right)
\]
columns, then with probability of at least $1-\delta$
\[
{\cal E}({\mathbf{x}}_{{\cal K},\mat R,k})\leq{\cal E}({\mathbf{x}}_{{\cal K},k})+\frac{(2\nu+\nu^{2})\TNormS{{\mathbf{f}}}}{n}
\]
where ${\mathbf{f}}$ is the expected value of ${\mathbf{b}}$ (recall the statistical
framework in section~\ref{subsec:statistical}).
\end{thm}

Before proving this theorem, we remark that the bound on the size
of the sketch is somewhat disappointing in the sense that it is useful
only if $d\gg n$ (since $\Trace{\ensuremath{{\bm{\mathrm{K}}}}}$ is likely to be large).
However, this is only a bound, and possibly a pessimistic one. Furthermore,
once the feature expanded data has been embedded in Euclidean space
(via \noun{TensorSketch}), it can be further compressed using standard
Euclidean space transforms like \noun{CountSketch} and subgaussian
maps (this is sometimes referred to as two-level sketching), or compression
can be applied from the left. We leave the task of improving the bound
and exploring additional compression techniques to future research.
\begin{proof}
The square singular values of $\Phi$ are exactly the eigenvalues
of $\ensuremath{{\bm{\mathrm{K}}}}$, so Theorem~\ref{thm:right_sketching} asserts that the
conclusions of the theorem hold if $\mat R^{\textsc{T}}$ provides $(\epsilon,\delta)$-approximate
Gram matrix for $\Phi$ where $\epsilon=O(\nu(\lambda_{k}-\lambda_{k+1})/\lambda_{1})$.
To that end, we combine the analysis of Avron et al.~\cite{ANW14}
of \noun{TensorSketch} with more recent results due to Cohen et al.~\cite{cohen2015optimal}.
Although not stated as a formal theorem, as part of a larger proof,
Avron et al. show that \noun{TensorSketch} has an OSE-moment property
that together with the results of Cohen et al.~\cite{cohen2015optimal}
imply that indeed the $(\epsilon,\delta)$-approximate Gram property
holds for the specified amounts of columns in $\mat R$.
\end{proof}

\section{\label{sec:experiments}Experiments}

In this section we report experimental results, on real data, that
illustrate and support the main results of the paper, and demonstrate
the ability of our algorithms to find appropriately regularized solutions.

\paragraph{Datasets. }

We experiment with three datasets, two regression datasets (\emph{Twitter
Buzz} and \emph{E2006-tfidf}) and one classification dataset (\emph{Gisette}).

\emph{Twitter Social Media Buzz~\cite{kawala2013predictions}} is
a regression dataset in which the goal is to predict the popularity
of topics as quantified by its mean number of active discussions given
77 predictor variables such as number of authors contributing to the
topic over time, average discussion lengths, number of interactions
between authors etc. We pre-process the data in a manner similar to
previous work \cite{lu2014fast,slawski2017compressed}. That is, several
of the original predictor variables, as well as the response variable
are log-transformed prior to analysis. We then center and scale to
unit norm. Finally, we add quadratic interactions which yielding a
total of $3080$ predictor variables (after preprocessing, the data
matrix is $583250\textrm{-by-}3080$). We used this dataset to explore
only sub-optimality of the objective and constraint satisfaction,
as we have found that the generalization error is very sensitive to
selection of the test set (when splitting a subset of the data to
training and testing).

\emph{E2006-tfidf}~\cite{kogan2009predicting} is regression dataset
where the features are extracted from SEC-mandated financial reports
published annually by a publicly traded company, and the quantity
to be predicted is volatility of stock returns, an empirical measure
of financial risk. We use the standard training-test split available
with the dataset\footnote{We downloaded the dataset from the LIBSVM website, \url{https://www.csie.ntu.edu.tw/~cjlin/libsvmtools/datasets/}.}.
We use this dataset only for testing generalization. The only pre-processing
we performed was subtracting the mean from the response variable,
and reintroducing it when issuing predictions.

The\emph{ Gisette} dataset is a binary classification dataset that
is a constructed from the MNIST dataset. The goal is to separate the
highly confusable digits '4' and '9'. The dataset has $6000$ data-points,
each having $5000$ features. We use the standard training-test split
available with the dataset (this dataset was downloaded from the same
website as the E2006-tfidf dataset). We convert the binary classification
problem to a regression problem using standard techniques (regularized
least squares classification). We use this dataset only for testing
generalization.

\paragraph{Baselines.}

A first reference are the performance of plain PCR. For small problems,
the dominant right subspace needed to compute the PCR solution can
be computed via MATLAB's dense SVD routine. For larger problems, we
compute the dominant right subspace using a PRIMME~\cite{SM10,WRS17},
a state-of-the-art iterative algorithm for SVD. As additional reference,
we also report results of two alternative algorithms: CLS and the
iterative algorithm of Frostig et al.~\cite{frostig2016principal}.
Both in the discussion, and in the graphs, we refer to the algorithm
Frostig et al.~ as ``Iterative-PCR''. We use the implementation
of Iterative-PCR supplied by the authors,\footnote{\url{https://github.com/cpmusco/fast-pcr}}
for which we used the default parameters, except for the ``tol''
parameter, which we set to $10^{-6}$ instead of the default $10^{-3}$.
We found that the use of tol=$10^{-3}$ produces results that generalize
poorly, while the use of tol=$10^{-6}$ produces much better results.
However, the running time of Iterative-PCR with tol=$10^{-6}$ is
considerably higher the the running time for tol=$10^{-3}$. Iterative-PCR
controls singular vector truncation via a cut-off parameter $\lambda$,
while in our experiments we set $k$ (the number of principal components
that are kept). We achieve this effect by setting $\lambda=(\sigma_{k}^{2}(\ensuremath{{\bm{\mathrm{A}}}})+\sigma_{k+1}^{2}(\ensuremath{{\bm{\mathrm{A}}}}))/2$
(when we report running times, we do not include the time to compute
the singular values). Finally, we remark that based of the documentation,
the algorithm analyzed by Frostig et al.~\cite{frostig2016principal}
does not completely correspond to the default parameters of the implementation
of Iterative-PCR supplied by the authors (e.g., the default parameter
for the ``method'' parameter is ``LANCZOS'', while ``EXPLICIT''
corresponds the algorithm analyzed in \cite{frostig2016principal}).

\paragraph{Sub-optimality of objective and constraint satisfaction.}

We explore the Twitter Buzz dataset from the optimization perspective,
namely measure the sub-optimality in the objective (vs. PCR) and constraint
satisfaction. Since $n\gg d$, we use left sketching with subgaussian
maps. We perform each experiment five times and report the median
value. Error bars, when present, represent the minimum and maximum
value of five runs. In the top panel, we use a fixed $k=60$ and vary
the sketch size (left sketching only), while in the bottom panel we
vary $k$ and set sketch size to be $s=4k$. The left panel explores
the value of the objective function, appropriately normalized (divided
by $\TNorm{\ensuremath{{\bm{\mathrm{A}}}}{\mathbf{x}}_{k}-{\mathbf{b}}}$ for fixed $k$, and divided by $\TNorm{{\mathbf{b}}}$
for varying $k$). The right panel explores the regularization effect
by examining the value of the constraints $\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\textsc{T}}{\mathbf{x}}_{k}}/\TNorm{{\mathbf{b}}}$.

In the left panel, we see that value of the objective for the sketched
PCR solution follows the value of objective for the PCR solution.
In general, as the sketch size increases, the variance in the objective
value reduces (top left graph). The normalized value of the constraint
for sketched PCR is rather small (as a reference we note that $\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\textsc{T}}{\mathbf{x}}_{OLS}}/\TNorm{{\mathbf{b}}}=0.4165$),
and generally decreases when the sketch size increases (top right
graph), but increases with $k$ for a fixed ratio between $s$ and
$k$ (bottom right graph). Furthermore, the results of sketched PCR
are very similar to the results of iterative PCR (bottom panel), while
running time is considerably shorter (see Table \ref{tab:time}).

The role of the constraints as a regularizer are illustrated by the
results for CLS (for fixed $k$ we use $t=4k$). As expected, CLS
achieves lower objective value at the price of larger constraint infeasibility.
The values of $\TNorm{\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}^{\textsc{T}}{\mathbf{x}}_{CLS}}/\TNorm{{\mathbf{b}}}$ are
much smaller than the OLS value, but much larger than the values for
sketched PCR. Furthermore, it is hard to control the regularization
effect for CLS: when sketch size increases the objective decreases
and the constraint increases (compare to PCR and sketched PCR, top
panel).

\begin{figure}
\noindent \begin{centering}
\begin{tabular}{ccc}
\includegraphics[width=0.45\textwidth]{twitter_objective_k_60} & ~ & \includegraphics[width=0.45\textwidth]{twitter_constraints_k_60}\tabularnewline
\includegraphics[width=0.45\textwidth]{twitter_objective_with_fpcr_slow} &  & \includegraphics[width=0.45\textwidth]{twitter_constraints_with_fpcr_slow}\tabularnewline
\end{tabular}
\par\end{centering}
\caption{\label{fig:twitter}Sub-optimality of objective and constraint satisfaction
for the Twitter Buzz dataset. }
\end{figure}

\paragraph{Generalization results.}

We also explored the prediction error and the tradeoffs between compression
and regularization. We perform each experiment five times and report
the median value. Error bars, when present, represent the minimum
and maximum value of those five runs.

We report the Mean Square Error (MSE) of predictions for the E2006.tfidf
dataset in Figure~\ref{fig:E2006}. We compare CLS, iterative-PCR,
right sketching and two-sided sketching (the matrix is too large for
exact PCR, and $d\gg n$ so right sketching is more appropriate).
In the left panel we fix $k=600$ and vary the sketch size. The MSE
decreases as the sketch size increases for both sketching methods.
For CLS, initially the MSE decreases and is close to the MSE of the
two sketching methods, but for large sketch sizes the MSE starts to
go up, likely due to decreased level of regularization. We note that
the minimum MSE achieved by CLS is larger than achieved by both sketching
methods. A similar phenomenon is observed when we vary the value of
$k$ in the right panel.

\begin{figure}
\noindent \begin{centering}
\begin{tabular}{ccc}
\includegraphics[width=0.45\textwidth]{E2006_prediction_fixed_k_with_pcr} & ~ & \includegraphics[width=0.45\textwidth]{E2006_prediction_varing_k_with_fpcr_slow}\tabularnewline
\end{tabular}
\par\end{centering}
\caption{\label{fig:E2006}Mean squared error of predictions for the E2006.tfidf
dataset. }
\end{figure}

We report the classification error for the Giesette dataset in Figure~\ref{fig:gisette}.
In the left panel we fix $k=400$ and vary the sketch size. As a reference,
the error rate of exact PCR (with $k=400$) is $2.8\%$ and the error
rate for OLS is 9.3\%. Left sketching has error rate very close to
the error rate of exact PCR, especially when $s$ is large enough.
Right sketching does not perform as well as left sketching, but it
too achieves low error rate for large $s$. For both methods, the
error rate drops as the sketch sizes increase, and the variance reduces.
For CLS the error rate and variance initially drops as the sketch
size increases, but eventually, when sketch size is large enough,
the error rate and the variance increases. This is hardly surprising:
as the sketch size increase, CLS approaches OLS. This is due to the
fact that CLS uses the compression to regularize, and when the sketch
size is large there is little regularization. In the right panel,
we vary the value of $k$ and set $s=4k$ (left sketching) and $t=4k$
(right sketching and CLS). Left sketch and PCR consistently achieve
about the same error rate. For small values of $t$, CLS performs
well, but when $t$ is too large the error starts to increase. In
contrast, right sketching continues to perform well with large values
of $k$. Again, we see that CLS mixes compression and regularization,
and one cannot use a large sketch size and modest amount of regularization
with CLS.

\begin{figure}
\noindent \begin{centering}
\begin{tabular}{ccc}
\includegraphics[width=0.45\textwidth]{gisette_classification_k_400} & ~ & \includegraphics[width=0.45\textwidth]{gisette_classification_varying_vs_with_fpcr_slow}\tabularnewline
\end{tabular}
\par\end{centering}
\caption{\label{fig:gisette}Classification error for the Gisette dataset. }
\end{figure}

\paragraph{Running time. }

In Table~\ref{tab:time} we report a sample of the various running
times of the different algorithms. All experiments were conducted
using MATLAB, although the sketching routines were written in C. Running
times were measured on a machine with a 6-core Intel Xeon Processor
E5-1650 v4 CPU and 128 GB of main memory, running Ubuntu 16.04. For
plain PCR, we report running time using PRIMME, which we ran with
default parameters and no preconditioner. For PRIMME, we cap the number
of iterations at 100,000, and write ``FAIL'' in the table if the
PRIMME failed to convergence within that cap. For iterative PCR we
also report running times when we set tol to the default value, and
reduce the max number of iteration from 40 to 10. This results in
much faster running time, but much degraded generalization (not reported),
e.g. for E2006 the test MSE for Iter-PCR (tol=$10^{-3}$, iter=10)
is 0.32. With respect to running time, iterative-PCR is competitive
with sketched PCR only for E2006, but with worse classification error.
Using PRIMME for PCR is not competitive with sketched PCR. However,
we stress that we experimented with only three datasets, so the comparison
is not comprehensive.

\begin{table}
\centering{}\caption{\label{tab:time}Running times (in seconds). For sketched PCR, we
report in brackets the type of sketching used (left, right, or two-sided).}
{\footnotesize{}}%
\begin{tabular}{|c|c|c|>{\centering}p{2.4cm}|>{\centering}p{2.4cm}|c|}
\hline
 & {\footnotesize{}CLS ($t=400)$} & {\footnotesize{}PRIMME-PCR} & {\footnotesize{}Iter-PCR }{\footnotesize \par}

{\footnotesize{}(tol=$10^{-3}$, iter=10)} & {\footnotesize{}Iter-PCR (tol=$10^{-6}$)} & {\footnotesize{}Sketched-PCR}\tabularnewline
\hline
\hline
{\footnotesize{}Twitter, $k=82$} & {\footnotesize{}21.2} & {\footnotesize{}1730} & {\footnotesize{}742} & {\footnotesize{}4067} & {\footnotesize{}4.7 (left)}\tabularnewline
\hline
{\footnotesize{}Twitter, $k=152$} & {\footnotesize{}48.9} & {\footnotesize{}5907} & {\footnotesize{}1278} & {\footnotesize{}7759} & {\footnotesize{}8.4 (left) }\tabularnewline
\hline
{\footnotesize{}E2006, $k=1000$} & {\footnotesize{}100} & {\footnotesize{}14694} & {\footnotesize{}140} & {\footnotesize{}601} & {\footnotesize{}150 (two sided)}\tabularnewline
\hline
{\footnotesize{}E2006, $k=2000$} & {\footnotesize{}270} & {\footnotesize{}FAIL} & {\footnotesize{}320} & {\footnotesize{}791} & {\footnotesize{}815 (two sided)}\tabularnewline
\hline
{\footnotesize{}Gisette, $k=400$} & {\footnotesize{}2.0} & {\footnotesize{}497} & {\footnotesize{}7.1} & {\footnotesize{}30.3} & {\footnotesize{}0.2 (left)}\tabularnewline
\hline
{\footnotesize{}Gisette, $k=1000$} & {\footnotesize{}9.3} & {\footnotesize{}FAIL} & {\footnotesize{}9.3} & {\footnotesize{}39.4} & {\footnotesize{}0.9 (left)}\tabularnewline
\hline
\end{tabular}
\end{table}

\section{Conclusions and Future work }

In this paper, we studied the use of sketching to accelerate the solution
of PCR and PCP. In particular, for a data matrix $\ensuremath{{\bm{\mathrm{A}}}}$, we relate
the PCR/PCP solution of $\ensuremath{{\bm{\mathrm{A}}}}\mat R$, where $\mat R$ is any dimensionality
reduction matrix, to the PCR/PCP solution of $\ensuremath{{\bm{\mathrm{A}}}}$. We presented
a notion of approximate PCR/PCP, motivated both from an optimization
perspective and from a statistical perspective, and provide conditions
on $\mat R$ that guarantee rigorous theoretical bounds. We then leverage
the aforementioned results to design fast, sketching based, algorithms
for approximate PCR/PCP, and demonstrate empirically the utility of
our proposed algorithms. Throughout, our focus in this paper has been
on algorithms that use the ``sketch-and-solve'' approach.

There are multiple ways in which the current work can be extended,
and the theoretical results improved. We have presented two notions
of approximation: approximate PCR and approximate PCP. Our results
for approximate PCR use only dimensionality reduction matrices $\mat R$
whose number of columns is equal to the target rank. It is natural
to conjecture that the use of dimensionality reduction matrices with
an higher number of columns will lead to stronger PCR bounds, but
we prove only PCP bounds. The underlying reason is that our bounds
for PCR are based on analyzing the distance between the column space
of $\mat R$ and the column space of $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$ . However,
once the number of columns in $\mat R$ is different from the number
of columns in $\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k}$, the definition of $d_{2}(\mat R,\ensuremath{{\bm{\mathrm{V}}}}_{\ensuremath{{\bm{\mathrm{A}}}},k})$
is no longer applicable. One possible strategy for analyzing PCR when
$\mat R$ has more than $k$ columns might be to use a generalization
of the distance between two subspaces that allows subspaces of different
size; see \cite{ye2016schubert} for such generalizations. Another
crucial component will then be to generalize the Davis-Kahan theorem
to bound such distances. We conjecture it is possible to derive algorithms
that depend on gaps between $\sigma_{k}$ and $\sigma_{k+l}$, where
$l$ is some oversampling parameter, as opposed to the smaller gap
between $\sigma_{k}$ and $\sigma_{k+1}$. We leave this for future
work.

Another interesting direction is in finding other ways to identify
a valid approximate dominant subspace. If we consider the statistical
perspective and inspect Eq.~\ref{eq:bias-var}, we see that all we
need is to find a subspace ${\cal S}\subseteq\range{\ensuremath{{\bm{\mathrm{A}}}}}$ of rank
$k$ such that $\FNorm{(\mat I-\ensuremath{{\bm{\mathrm{P}}}}_{{\cal S}})\ensuremath{{\bm{\mathrm{A}}}}}$ is small,
while our theoretical results try to achieve a stronger bound: having
the dominant subspaces align. One possible way for finding such a
${\cal S}$ is using so-called Projection-cost Preserving Sketches~\cite{CohenEMMP15}.
We leave this for future work.

\subsection*{Acknowledgments.}

This research was supported by the Israel Science Foundation (grant
no. 1272/17) and by an IBM Faculty Award.

\bibliographystyle{plain}

redpajama_set_name: RedPajamaArXiv