\section*{Acknowledgment}
A. Zouzias would like to thank Mark Rudelson for several discussions on the approximate matrix multiplication problem. This research has received funding from the ERC under the EU's Seventh Framework Programme (FP7/2007-2013) / ERC grant agreement $n^{o}$ 259569.


\bibliographystyle{IEEEtran}

\section{Application to data-driven low-dimensional embedding}
As an application of the result above, we consider the following question:

\vskip.1in
\noindent \textsc{Problem:} {\it 
Given a collection $n$ points in $\mathbb{R}^d$, forming a matrix $\mat{A} \in \mathbb{R}^{n \times d}$ and an error parameter $\varepsilon>0$, construct efficient and approximately accurate low-dimensional embedding $\mat{G} \in \mathbb{R}^{t \times d}$ such that:}
\begin{align}
\|\mat{G}(\mathbf{a}_i - \mathbf{a}_j)^T\|_2^2 \stackrel{\pm \varepsilon}{\approx} \|\mathbf{a}_i - \mathbf{a}_j\|_2^2,\quad \forall \mathbf{a}_i \in \mathbb{R}^{1 \times d}, ~\mathbf{a}_i \in \text{\texttt{Rows}}(\mat{A}). \nonumber
\end{align}
\noindent {\it The target dimension $t$ is to heavily depend on the input data matrix $\mat{A}$ and be independent of the input dimensions.}
\vskip.1in


Inspired by Theorem~\ref{thm:main}, the following theorem proposes a data-dependent randomized low-dimensional embedding $\mat{G}$ with the following guarantees:
\begin{theorem}\label{eq:main_thm}
Let $\mat{A} \in \mathbb{R}^{n \times d}$. Moreover, assume that $\mat{A}$ with nuclear rank $ \nr{\mat{A}} \ll O(\log n)$. If $t = \Omega\left( \frac{\nr{\mat{A}} + \log( \log(1/\varepsilon)) + \log(1/\delta)}{\varepsilon^2}\right)$, then the following inequalities hold 
\begin{align}
\forall i, j \in [n]\quad \left| \left\|\mat{G}\left(\mathbf{a}_i - \mathbf{a}_j\right)^\top\right\|_2^2 - \left\|\mathbf{a}_i - \mathbf{a}_j\right\|_2^2 \right| \leq 2\varepsilon \norm{\mat{A}}^2, \nonumber
\end{align}
with probability at least $1-\delta$.
\end{theorem} 

To fully specify the algorithmic procedure followed in practice, only $\mat{A} \in \mathbb{R}^{n \times d}$ and constants $\delta,~\varepsilon \in (0, 1)$ are given as input and $\mat{G} \in \mathbb{R}^{t \times d}$ is returned. Moreover, $\mat{G}$ guarantees to preserve the distances between the rows of $\mat{A}$ with additive error $2\varepsilon \norm{\mat{A}}^2$ and with probability at least $1- \delta$. The exact steps followed are given in Algorithm \ref{algo:1}.

\begin{algorithm}[!htp]
\caption{Data-Driven Random Linear Embedding}\label{algo:1}
\begin{small}
\begin{algorithmic}[1]
\Statex ~~~\textbf{Input:} $\mat{A} \in \mathbb{R}^{n \times d}$, $\varepsilon, \delta \in (0, 1)$
\algrule
\State ~~~ Compute (or approximate) $\|\mat{A}\|_{\star}$ and $\norm{\mat{A}}$.
\State ~~~ Set parameter $t = \Omega\left( \frac{\nr{\mat{A}} + \log \log(1/\varepsilon) + \log(1/\delta)}{\varepsilon^2} \right)$.
\State ~~~ Generate $\mat{G} \in \mathbb{R}^{t \times n}$ according to Theorem \ref{eq:main_thm}.
\algrule
\State ~~~\textbf{Output:} Matrix $\mat{G}$ and $\mat{G}\mat{A}$ that satisfies Theorem \ref{eq:main_thm} with probability at least $1 - \delta$.
\end{algorithmic}
\end{small}
\end{algorithm}

\begin{proof}
By substituting $\mat{X} = \mat{Y}^\top = \mat{A}$ in Theorem \ref{thm:main}, we have:
\begin{align}
\mathbb{P}\left( \norm{\mat{A} \mat{G}^\top \mat{G} \mat{A}^\top - \mat{A} \mat{A}^\top} \leq \varepsilon \norm{\mat{A}}^2\right) \geq 1 - \delta. \nonumber
\end{align}
By the definition of the spectral norm, the condition above can be further written as the following maximization problem:
\begin{align}
\max_{\mathbf{x}: \mathbf{x} \in \mathbb{R}^n, \|\mathbf{x}\|_2 = 1} \left\{ \left| \mathbf{x}^\top \left( \mat{A} \mat{G}^\top \mat{G} \mat{A}^\top - \mat{A} \mat{A}^\top\right) \mathbf{x} \right| \right\} \leq \varepsilon \norm{\mat{A}}^2. \nonumber
\end{align} Moreover, it is obvious that, if we restrict the search space $\mathbf{x} \in \mathbb{R}^n$ to $\mathbf{x} \in \mathbb{B}^n$ where:
\begin{align}
\mathbb{B}^n = \left\{ \mathbf{y} ~ | ~\mathbf{y} = \frac{1}{\sqrt{2}}\left(\mathbf{e}_i - \mathbf{e}_j\right), \forall i, j \in [n] \right\}, \nonumber
\end{align} where $\mathbf{e}_i$ is the standard basis vector with 1 in the $i$-th position, we further have:
\begin{align}
\max_{\mathbf{x} : \mathbf{x} \in \mathbb{B}^n} &\left\{ \left| \mathbf{x}^\top \left( \mat{A} \mat{G}^\top \mat{G} \mat{A}^\top - \mat{A} \mat{A}^\top\right) \mathbf{x} \right| \right\} \nonumber \\ &\leq \max_{\mathbf{x}: \mathbf{x} \in \mathbb{R}^n, \|\mathbf{x}\|_2 = 1} \left\{ \left| \mathbf{x}^\top \left( \mat{A} \mat{G}^\top \mat{G} \mat{A}^\top - \mat{A} \mat{A}^\top\right) \mathbf{x} \right| \right\} \nonumber
\end{align} which leads to:
\begin{align}{\label{eq:2}}
\mathbb{P}\left(\max_{\mathbf{x} : \mathbf{x} \in \mathbb{B}^n} \left\{ \left| \mathbf{x}^\top \left( \mat{A} \mat{G}^\top \mat{G} \mat{A}^\top - \mat{A} \mat{A}^\top\right) \mathbf{x} \right| \right\} \leq \varepsilon \norm{\mat{A}}^2\right) \geq 1 - \delta.
\end{align} 
Observe also that $\left| \mathbf{x}^\top \left( \mat{A} \mat{G}^\top \mat{G} \mat{A}^\top - \mat{A} \mat{A}^\top\right) \mathbf{x} \right| = \left| \|\mat{G} \mat{A}^\top\mathbf{x}\|_2^2 - \|\mat{A}^\top \mathbf{x}\|_2^2 \right|$, which further transforms \eqref{eq:2} as:
\begin{align}{\label{eq:3}}
\mathbb{P}\left(\max_{\mathbf{x} : \mathbf{x} \in \mathbb{B}^n} \left\{ \left| \|\mat{G} \mat{A}^\top\mathbf{x}\|_2^2 - \|\mat{A}^\top \mathbf{x}\|_2^2 \right| \right\} \leq \varepsilon \|\mat{A}\|_2^2\right) \geq 1 - \delta.
\end{align} For any vector $\mathbf{x} \in \mathbb{B}^n$, we observe that $\mathbf{x}^\top \mat{A} \mat{G}^\top = 1/\sqrt{2}\left(\mathbf{a}_i - \mathbf{a}_j\right)\mat{G}^\top$ while $\mat{G} \mat{A}^\top \mathbf{x} = 1/\sqrt{2}\mat{G} \left(\mathbf{a}_i - \mathbf{a}_j\right)^\top$. Similarly, $\mathbf{x}^\top \mat{A}= 1/\sqrt{2}\left(\mathbf{a}_i - \mathbf{a}_j\right)$ and $\mat{A}^\top \mathbf{x} = 1/\sqrt{2}\left(\mathbf{a}_i - \mathbf{a}_j\right)^\top$. Thus, \eqref{eq:3} becomes:
\begin{small}
\begin{align}{\label{eq:4}}
\mathbb{P}\Big(\max_{i, j \in [d], i \neq j} \left\{ \left|~\left\|\mat{G}\left(\mathbf{a}_i - \mathbf{a}_j\right)^\top\right\|_2^2 - \left\|\mathbf{a}_i - \mathbf{a}_j\right\|_2^2~\right| \right\} &\leq 2 \varepsilon \|\mat{A}\|_2^2\Big) \nonumber \\ &\geq 1 - \delta.
\end{align}
\end{small} Since \eqref{eq:4} is satisfied for the maximizing combination of $i, j$, we can safely remove the maximization to get:
\begin{align}
\mathbb{P}\left(\left|~\left\|\mat{G}\left(\mathbf{a}_i - \mathbf{a}_j\right)^\top\right\|_2^2 - \left\|\mathbf{a}_i - \mathbf{a}_j\right\|_2^2~\right| \leq 2\varepsilon \|\mat{A}\|_2^2\right) \geq 1 - \delta, \nonumber
\end{align} which completes the proof.
\end{proof} A complete set of experiments will be included in an extended version of the paper.
\section*{Appendix}

In our analysis, we make use of the following concentration bound on the operator norm of the product of a fixed matrix with a Gaussian matrix which is a direct consequence of concentration of Lipschitz function on Gaussian space, see e.g.~\cite[p.~10]{zouzias2013randomized}.
\begin{lemma}\label{lem:normGaussian}
Let $\mat{X} \in \mathbb{R}^{n\times d}$ and $t\geq 1$. For every $\tau>0$:
\begin{equation}
\Prob{ \norm{\mat{X}\mat{G}^\top} \geq \frobnorm{\mat{X}} / \sqrt{t} + \norm{\mat{X}} + \tau \norm{\mat{X}}/\sqrt{t}} \leq \exp(-\tau^2 /8). \nonumber
\end{equation}
\end{lemma}


\subsection{Proof of Lemma~\ref{lem:prob_bound}}
In all the cases below, observe $\gamma_{l,s} > \varepsilon$ by definition. Using Eqn.~\eqref{eq:proba3}, we consider the following two cases:
\begin{itemize}
\item [$(i)$] $\gamma_{l,s}  < 1$: in this case, we have $\min(\gamma_{l,s}, 1) \geq \min(\varepsilon, 1) = \varepsilon$. Thus, \eqref{eq:proba3} becomes:
\begin{align}
\mathbb{P}\left(\Omega(l, s)\right) &\leq c_2^{r(\mat{X}_{\theta}^{l}) + r(\mat{Y}_{\theta}^s)}\cdot e^{-c_1\left(\varepsilon + c_2\frac{r(\mat{X}^{l}_{\theta}) + r(\mat{Y}^{s}_{\theta})}{ c_1 \varepsilon t}\right)t \cdot \varepsilon} \nonumber \\ 
&\leq e^{-c_1\varepsilon^2 t} \cdot e^{-c_2\left(r(\mat{X}^{l}_{\theta}) + r(\mat{Y}^{s}_{\theta})\right)} \cdot e^{\ln (c_2)\left(r(\mat{X}^{l}_{\theta}) + r(\mat{Y}^{s}_{\theta})\right)} \nonumber \\
&\leq e^{-c_1\varepsilon^2 t}
\end{align}
\item [$(ii)$] $\gamma_{l,s}  \geq  1$: in this case, $\min(\gamma_{l,s}, 1) = 1$. Following the same steps as above, we have:
\begin{align}
\mathbb{P}\left(\Omega(l, s)\right) &\leq c_2^{r(\mat{X}_{\theta}^{l}) + r(\mat{Y}_{\theta}^s)}\cdot e^{-c_1\left(\varepsilon + c_2\frac{r(\mat{X}^{l}_{\theta}) + r(\mat{Y}^{s}_{\theta})}{ c_1 \varepsilon t}\right)t} \nonumber \\ 
&\leq e^{-c_1\varepsilon t} \cdot e^{\left(1 - \frac{1}{\varepsilon}\right)c_2\left(r(\mat{X}^{l}_{\theta}) + r(\mat{Y}^{s}_{\theta})\right)}.  \nonumber
\end{align}
However, since $0 < \varepsilon < 1$, we have $1 - \frac{1}{\varepsilon} < 0$ and thus, the above inequality can be further upper bounded by:
\begin{align}
\mathbb{P}\left(\Omega(l, s)\right) \leq e^{-c_1\varepsilon t}. \nonumber
\end{align}
\end{itemize}
\subsection{Proof of Lemma~\ref{lem:tailTerms}}
Apply Lemma~\ref{lem:normGaussian} with $\tau = \sqrt{t}$ on all four matrices $\mat{X}_{\theta}^c$, $\mat{X}_{\theta}$, $\mat{Y}_{\theta}$ and $\mat{Y}_{\theta}^c$. With probability at least $1-\exp(-t/8)$:
\begin{align}
		\norm{\mat{X}_{\theta}^c\mat{G}^\top} &\leq \frac{\frobnorm{\mat{X}_{\theta}^c}}{\sqrt{t}}  + \left(1 + \frac{\tau}{\sqrt{t}}\right)\norm{\mat{X}_{\theta}^c} \nonumber \\ &\stackrel{\tau = \sqrt{t}}{=} \frac{1}{\sqrt{t}} \frobnorm{\mat{X}_{\theta}^c} + 2\varepsilon \nonumber \\
										&\leq \frac{1}{\sqrt{t}} \frobnorm{\mat{X}} + 2\varepsilon \leq \frac{1}{\sqrt{t}} \cdot \|\mat{X}\|_{\star} + 2\varepsilon  \nonumber\\
									   &\stackrel{t \geq \theta}{\leq} \frac{1}{\sqrt{\theta}} \nr{\mat{X}} + 2\varepsilon \leq \frac{\varepsilon}{\sqrt{c_3}} + 2\varepsilon \stackrel{c_3 > 1}{\leq} 3\varepsilon \nonumber
\end{align} 
by definition of $\theta$ and $\norm{\mat{X}_{\theta}^c} \leq \varepsilon$.
Similarly, with probability at least $1-\exp(-t/8)$:
\begin{align}
\norm{\mat{X}_{\theta}\mat{G}^\top} &\leq  \frac{\frobnorm{\mat{X}_{\theta}}}{\sqrt{t}} + \left(1 + \frac{\tau}{\sqrt{t}}\right)\norm{\mat{X}_{\theta}} \nonumber \\ &\stackrel{\tau = \sqrt{t}}{=} \frac{1}{\sqrt{t}} \frobnorm{\mat{X}_{\theta}} + 2 \nonumber \\
								&\leq \frac{1}{\sqrt{t}} \frobnorm{\mat{X}} + 2 \leq \frac{1}{\sqrt{t}} \cdot \|\mat{X}\|_{\star} + 2  \nonumber\\
	     					   &\stackrel{t \geq \theta}{\leq}  \frac{1}{\sqrt{\theta}} \nr{\mat{X}} + 2\nonumber \leq \frac{\varepsilon}{\sqrt{c_3}} + 2 \stackrel{c_3 > 1}{\leq} 2(\varepsilon + 1) \nonumber
\end{align}
		by definition of $\theta$ and $\norm{\mat{X}_{\theta}} = 1$. 
We further observe $\norm{\mat{Y}_{\theta}^c\mat{G}^\top} \leq 3\varepsilon$ and $\norm{\mat{X}_{\theta}\mat{G}^\top} \leq 2(\varepsilon + 1)$ both with probability at least $1-\exp(-t/8)$. 
	Union bound all the above four results, it follows that
	\begin{align}
	&\norm{ \mat{X}_{\theta}^c\mat{G}^\top \mat{G}\mat{Y}_{\theta}} + \norm{\mat{X}_{\theta}\mat{G}^\top \mat{G}\mat{Y}_{\theta}^c} + \norm{ \mat{X}_{\theta}^c\mat{G}^\top \mat{G}\mat{Y}_{\theta}^c}  \nonumber \\ 
	&\leq 12\varepsilon\left(\varepsilon + 1\right) + 9\varepsilon^2 \leq 12\varepsilon  + 21\varepsilon^2 \leq 33\varepsilon \nonumber
	\end{align} 
which holds with probability at most $4\cdot e^{-t/8}$.  By definition of $t$, the above inequality is violated with probability at most $\delta/2$.	

\section{Conclusions}

We present a novel analysis for a class of randomized and provably $\varepsilon$-accurate algorithms for the problem of matrix multiplication. 
As an application of this result, we show the utilization of the proposed scheme on data-driven low dimensional embeddings with additive error approximation, in the case where the data live on a subspace characterized by a small \emph{nuclear rank}. 

An interesting question to pursue lies in the substitution of nuclear rank by stable rank: recent developments on this topic \cite{zouzias2013randomized} show similar results using the latter metric as the intrinsic data dimension; a weaker assumption than the nuclear rank. However, the dependence on the approximation error is of the order $O(\frac{1}{\varepsilon^4})$, as opposed to $O(\frac{1}{\varepsilon^2})$ presented in this work. We hope this paper triggers future efforts to improve stable rank-based bounds with respect to error dependency.
\section{Our Contributions}
Here, we highlight the main contributions of this manuscript:
\begin{enumerate}
\item [$(i)$] We provide a novel analysis on the approximate matrix multiplication problem based on the notion of nuclear rank (See Theorem~\ref{thm:main}).
\item [$(ii)$] We demonstrate an application of the approximate matrix multiplication bound to dimensionality reduction with relative guarantees.
\end{enumerate}
To highlight the potential of the second contribution, consider the scenario of an almost one-dimensional point set, i.e., randomly generated points centred around a line with small variance in high dimensional space. In this setting, intuitively we expect that a constant number of dimensions are necessary in order to preserve the geometry of the dataset. However, the JL lemma~\cite{johnson1984extensions} (pessimistically) proposes that logarithmic number of dimensions is necessary, which contradicts out expectations for less dimensions.



\section{Experiments}
\section{Introduction}

Living in the era of \emph{Big Data}, the excess of available information constitutes its manipulation and interpretation a strenuous task: In contrast to conventional wisdom where more data is a source of ``simplicity'' in statistical terms \cite{chandrasekaran2013computational}, large datasets embody preprocessing tasks of high time- and space-complexity, jeopardizing any hope for data analysis within reasonable time and with low computational cost. Due to such difficulties, one might be interested in \emph{screening} the data, even if \emph{no prior information is available}. I.e., identifying a coreset such that most of the latent structure/information that we want to infer is maintaned within an error that we can control. Such approaches have been witnessed in a broad class of data analysis problems such as data matrix sparsification for accelerated spectral calculation \cite{achlioptas2013near}, feature selection (a.k.a. column subset selection) for better interpretation of the results \cite{mahoney2009cur} and low-complexity calculations \cite{BoutsidisMD09}, feature extraction via dimensionality reduction techniques \cite{indyk1998approximate, BoutsidisZD10}, etc. 


In this work, we analyze a particular approximation algorithm of~\cite{sarlos2006improved} for the task of matrix multiplication with respect to their spectral norm \cite{zouzias2013randomized}. Within this context, we focus on approximately computing the product of two matrices \emph{when their intrinsic dimensionality is low, as expressed by the nuclear rank.} Here, we define the nuclear rank as the ratio of the nuclear norm of a matrix over its spectral norm. We provide an elementary proof based on the randomized algorithm described in \cite{sarlos2006improved,zouzias2013randomized}. As a result, we further strengthen the performance of the proposed scheme in situations where the allowed approximation error is small: Using a weaker notion of intrinsic dimensionality (i.e., nuclear rank instead of stable rank; see the definitions later in text), the dependence on the approximation error $\varepsilon$ is improved to $\mathcal{O}(1/\varepsilon^2)$ instead of $\mathcal{O}(1/\varepsilon^4)$, indicating that the proposed scheme scales better when $\varepsilon$ decreases. Table~\ref{table:mm} places our result into context with prior works.
\begin{table}[!h]
\centering
\caption{Summary of results on $\varepsilon$-approximate matrix multiplication $\mat{A} \cdot \mat{B}$ with respect to spectral norm -- $\nr{\cdot}$ and $\sr{\cdot}$ denote the nuclear and stable rank of a matrix, respectively.}\label{table:mm}
\ra{1.3}
\begin{tabular}{l c c c c} \toprule
\textbf{Metric} & \phantom{ab} & \textbf{\# of dimensions} & \phantom{ab} & \textbf{Reference} \\
\cmidrule{1-1} \cmidrule{3-3} \cmidrule{5-5}
 Rank		& & $\mathcal{O}( ( r(\mat{A}) + r(\mat{B}) ) / \varepsilon^2)$	& & \cite{MZ11}	\\
 Stable rank	& & $\mathcal{O}( (\sr{\mat{A}} + \sr{\mat{B}}) / \varepsilon^4)$	& 	& \cite{zouzias2013randomized}	\\
 Nuclear norm	&  & $\mathcal{O}( (\nr{\mat{A}} + \nr{\mat{B}}) / \varepsilon^2)$   &	&   Theorem~\ref{thm:main} \\
\bottomrule
\end{tabular}
\end{table}


As an application, we use this result to design dimensionality reducing linear embeddings that operate on a given dataset and preserve Euclidean distances between data points. In general, it is well-known that simply rotating, scaling and translating in a random way is adequate for this task: According to JL lemma \cite{johnson1984extensions}, such linear mapping approximately preserves (i.e., within a $\varepsilon$-radius of points) the distances between all the pairs of points. While there are many random constructions that achieve this property with high probability, such schemes are oblivious to the geometry of the data set at hand, leaving space for further improvements.

\vskip.07in
\noindent \textbf{Contributions:} The main contributions of this manuscript are the following:
\begin{enumerate}
\item [$(i)$] We provide a novel analysis on the approximate matrix multiplication problem based on the notion of nuclear rank (See Theorem~\ref{thm:main}).
\item [$(ii)$] We demonstrate an application of the approximate matrix multiplication bound to dimensionality reduction with relative guarantees.
\end{enumerate}




\section{Preliminaries}

A scalar is denoted by an italic letter, e.g. $a$. A column vector is denoted by a bold lowercase letter, e.g. $\mathbf{a} \in \mathbb{R}^n$ whose $i$-th entry is $a_i$. A matrix is denoted by a mathtype uppercase letter, e.g., $\mat{A} \in \mathbb{R}^{n \times d}$ with $(i, j)$-th entry $\mat{A}_{i, j}$. To denote the $i$-th row and $j$-th column of $\mat{A}$, we use the $\mathbf{a}_i$ and $\mathbf{a}^{j}$, respectively. Use $\norm{\mat{A}} = \max_{\mathbf{x}: \|\mathbf{x}\|_2 = 1} \|\mat{A}\mathbf{x}\|_2$ to denote the spectral norm of $\mat{A}$ and $\frobnorm{\mat{A}} = \sqrt{\sum_i \sum_j \mat{A}_{i, j}^2}$ to represent its Frobenius norm. 

In our analysis, we utilize the notion of {\it nuclear norm}. Given a matrix $\mat{A} \in \mathbb{R}^{n \times d}$ with rank $r(\mat{A})$, its nuclear norm $\|\mat{A}\|_{\star}$ is given by $\|\mat{A}\|_{\star} = \sum_{i = 1}^{r(\mat{A})} \sigma_i$ where $\sigma_i$ is the $i$-th largest singular value; for multiple matrices, we also use $\xi_i$ to denote a singular value. Using the nuclear norm, we define the \emph{nuclear rank} of $\mat{A}$ as $\nr{\mat{A}} := \|\mat{A}\|_{\star} / \norm{\mat{A}}$. Moreover, we use the notion of \emph{stable rank}, which is defined as $\sr{\mat{A}} = \frobnorm{\mat{A}}^2 / \norm{\mat{A}}^2$.

Throughout the paper, we will denote by $\mat{G}$ a $t \times d$ random matrix whose entries are independent Gaussian random variables with variance $1/t$, i.e., $\mat{G}_{ij} \sim \mathcal{N}(0, 1/t)$.

\section{Related Work}
There is rich literature on the topic of linear dimensionality reduction and the techniques utilized for this purpose, mainly due to the wide range of applications it covers. Here, we highlight a few approaches that provide bounds depending on the input data: These bounds are usually stronger than the classic JL lemma, provided that the input points have low ``geometric complexity'', i.e., points that lie on a low dimensional subspace~\cite{sarlos2006improved}, lie on a manifold~\cite{baraniuk,clarkson2008tighter}, etc. Here, we focus only on results that provide theoretical guarantees.

For the special case where the input points lie on the unit sphere, there is a close connection between the Talagrand functional, denoted as $\gamma_2$, and the required number of dimensions for point distance preservation~\cite{mendelson},\cite{Klartag}. These results can be viewed as stronger bounds compared to the JL lemma. Baraniuk et al. provide bounds on the number of dimensions required by a random linear embedding to preserve Euclidean distances for the case of smooth manifolds~\cite{baraniuk}; for improvements and a nice exposition on the topic, see~\cite{clarkson2008tighter}. Sarlos provided improved bounds for algorithms on large matrices with applications on matrix multiplications, linear regression and low rank matrix approximation~\cite{sarlos2006improved}. On a similar flavor, the authors in~\cite{related:movingPoints} provide bounds for relative error approximations of points that lie on a surface. Several approaches exist for non-linear dimensionality reduction as well~\cite{JL:beyond,JL:beyond2}. From a different perspective, Indyk and Naor consider the ``doubling dimension'' and other related measures of dimensionality for preserving nearest neighbor queries~\cite{nnEmbed}.

Recent developments in \cite{hegde2012convex} describe \emph{deterministic constructions} of linear embeddings in polynomial time via SemiDefinite Programming (SDP) relaxations, under the assumption that the data is known \emph{apriori} and fixed. 
For the latest developments on this topic, we refer the reader to \cite{grant2013nearly}.


\section{Main result}
We consider the following problem:
\vskip.1in
\noindent \textsc{Approximate Matrix Multiplication Problem:} Let $\mat{X} \in \mathbb{R}^{n \times d}$ and $\mat{Y} \in \mathbb{R}^{d \times m}$ be two arbitrary matrices, $\varepsilon > 0$ is an approximation parameter and $0<\delta<1$ the failure probability. We desire to construct \emph{sketched} matrices $\widehat{\mat{X}} \in \mathbb{R}^{n \times t}$ and $\widehat{\mat{Y}} \in \mathbb{R}^{t \times m}$ where $t \ll d$ such that:
\begin{align}\nonumber
\norm{\widehat{\mat{X}} \widehat{\mat{Y}} - \mat{X} \mat{Y}} \leq  \varepsilon \norm{\mat{X}} \cdot \norm{\mat{Y}},
\end{align} holds with probability at least $1-\delta$.
\vskip.1in
We now state the main theorem of the paper:
\begin{theorem}{\label{thm:main}}
Fix $0 < \varepsilon, \delta < 1$ and assume arbitrary matrices $\mat{X} \in \mathbb{R}^{n \times d}$ and $\mat{Y} \in \mathbb{R}^{d \times m}$. Set $\widehat{\mat{X}} = \mat{X} \mat{G}^\top$ and $\widehat{\mat{Y}} = \mat{G} \mat{Y}$. If $t = \Omega \left( \frac{\nr{\mat{X}} + \nr{\mat{Y}} + \log(\log(1/\varepsilon)) +\log(1/\delta)} {\varepsilon^2} \right)$, then the following holds:
\begin{equation}{\label{eq:thm}}
\Prob{ \norm{\widehat{\mat{X}} \widehat{\mat{Y}} - \mat{X} \mat{Y}} \leq  \varepsilon \norm{\mat{X}} \cdot \norm{\mat{Y}}} \geq 1 - \delta,
\end{equation}
\end{theorem}
We devote the rest of this section to prove Theorem~\ref{thm:main}. We recall a well-established result from the literature for Gaussian matrices; observe the lack of any upper bound on the error parameter.
\begin{lemma}\label{lem:mm:rank}
Fix $\zeta > 0$ and let $\mat{X} \in \mathbb{R}^{n \times d}$ and $\mat{Y} \in \mathbb{R}^{d \times m}$. Set $\widehat{\mat{X}} = \mat{X} \mat{G}^\top$ and $\widehat{\mat{Y}} = \mat{G} \mat{Y}$. Then,
\begin{small}
\begin{align}\nonumber
\Prob{ \norm{\widehat{\mat{X}}\widehat{\mat{Y}} - \mat{X} \mat{Y} } > \zeta \norm{\mat{X}}\norm{\mat{Y}}} \leq c_2^{r(\mat{X})+r(\mat{Y})} e^{- c_1t \zeta \cdot \min(\zeta, 1)}
\end{align}
\end{small}
where $c_1$ is the Hanson-Wright constant~\cite{ineq:HR} and $c_2=18$.
\end{lemma}
\begin{proof}
The proof is a corollary of the Hanson-Wright inequality~\cite[Theorem 1.1]{ineq:HR}, combined with a dense net argument on the unit sphere defined by the union of the column and row span of $\mat{Y}$ and $\mat{X}$, respectively.
\end{proof}

By homogeneity, we observe that, to prove $\norm{\mat{X} \mat{G}^\top \mat{G} \mat{Y} - \mat{X} \mat{Y}} \leq \varepsilon \norm{\mat{X}} \norm{\mat{Y}}$ is satisfied with some probability, it suffices to prove the same argument for $\norm{\frac{\mat{X}}{\norm{\mat{X}}} \mat{G}^\top \mat{G} \frac{\mat{Y}}{\norm{\mat{Y}}} - \frac{\mat{X}}{\norm{\mat{X}}} \cdot \frac{\mat{Y}}{\norm{\mat{Y}}}} \leq \varepsilon$. Thus, without loss of generality, we can assume that $\norm{\mat{X}} = \norm{\mat{Y}} = 1$.

Let $\mat{X} = \sum_{i=1}^{r(\mat{X})} \sigma_i \mathbf{u}_i \mathbf{v}_i^\top$ be the singular value decomposition (SVD) of $\mat{X}$, where $\sigma_1, \dots, \sigma_{r(\mat{X})}$ denote the singular values of $\mat{X}$ and, $\mathbf{u}_1, \dots, \mathbf{u}_{r(\mat{X})} \in \mathbb{R}^n$ and $\mathbf{v}_1, \dots, \mathbf{v}_{r(\mat{X})} \in \mathbb{R}^d$ denote the left and right singular vectors, respectively. Similarly, we define the SVD of $\mat{Y}$ as $\mat{Y} = \sum_{j=1}^{r(\mat{Y})} \xi_j \mathbf{p}_j \mathbf{q}_j^\top$.

Define $\theta := \left \lfloor c_3\cdot \frac{\nr{\mat{X}} + \nr{\mat{Y}}}{\varepsilon^2}\right\rfloor$ where $c_3 > 1$ is a constant and set $t:= \theta + 8\ln(8/\delta) /\varepsilon^2 + \ln\left( \left\lceil \ln(e/\varepsilon)\right\rceil\right) / (c_1 \cdot \varepsilon^2)$. Given $\theta$, one can decompose $\mat{X}$ as $\mat{X} = \mat{X}_{\theta} + \mat{X}_{\theta}^c$ where $\mat{X}_{\theta}$ represents the \emph{best} rank-$\theta$ approximation of $\mat{X}$ and $\mat{X}_{\theta}^c := \mat{X} - \mat{X}_{\theta}$. Similarly, we can decompose $\mat{Y} = \mat{Y}_{\theta} + \mat{Y}_{\theta}^c$.

Now, by the triangle inequality:
\begin{align}
\|\mat{X} \mat{G}^\top \mat{G} \mat{Y} &- \mat{X} \mat{Y}\|_2 \leq \norm{ \mat{X}_{\theta}\mat{G}^\top \mat{G}\mat{Y}_{\theta} - \mat{X}_{\theta} \mat{Y}_{\theta}} \label{eq:result_1}\\
&+ \norm{\mat{X}_{\theta}^c\mat{G}^\top \mat{G}\mat{Y}_{\theta}} + \norm{\mat{X}_{\theta}\mat{G}^\top \mat{G}\mat{Y}_{\theta}^c} + \norm{\mat{X}_{\theta}^c\mat{G}^\top \mat{G}\mat{Y}_{\theta}^c} \label{eq:termNormGaussian}\\
&+ \norm{\mat{X}_{\theta}^c\mat{Y}_{\theta}} + \norm{\mat{X}_{\theta}\mat{Y}_{\theta}^c} + \norm{\mat{X}_{\theta}^c\mat{Y}_{\theta}^c}. \label{eq:termConstant}
\end{align} 

Bounding the term appearing on the right hand side of \eqref{eq:result_1} is the most challenging task and it is the main technical contribution of this paper. Before we start, we need to define some notation. Let $l, s$ be non-negative integers and define the sets $\mathcal{I} = \lbrace \mathcal{I}_1, \mathcal{I}_2, \dots \rbrace$ and $\mathcal{J} = \lbrace \mathcal{J}_1, \mathcal{J}_2, \dots \rbrace$ where:
\begin{align}{\label{eq:def_sets1}} \vspace{-0.1cm}
\mathcal{I}_l := \left\{ i ~ : ~e^{-l} \leq \sigma_i < e^{-l + 1} \right\}, \vspace{-0.1cm}
\end{align} and
\begin{align}{\label{eq:def_sets2}} \vspace{-0.1cm}
\mathcal{J}_s := \left\{ j ~ : ~e^{-s} \leq \xi_j < e^{-s + 1} \right\}. \vspace{-0.1cm}
\end{align}
We remind that, since $\norm{\mat{X}} = \norm{\mat{Y}} = 1$, $\sigma_i \leq 1, ~\xi_j \leq 1$, for any $i, j$.

Given a set $\mathcal{I}_l$, we define $\mat{X}^{l}_{\theta} \in \mathbb{R}^{n \times d}$ as $\mat{X}^{l}_{\theta} := \sum_{i \in \mathcal{I}_l} \sigma_i \mathbf{u}_i \mathbf{v}_i^\top$; similarly, we have $\mat{Y}^{s}_{\theta} := \sum_{j \in \mathcal{J}_s} \xi_j \mathbf{p}_j \mathbf{q}_j^\top$ for the case of $\mat{Y}$. We highlight the following key observations:
\begin{itemize}
\item [$(i)$] $\mat{X}_{\theta} = \sum_{l = 1}^{L} \mat{X}^{l}_{\theta}$ and $\mat{Y}_{\theta} = \sum_{s = 1}^{S} \mat{Y}^{s}_{\theta}$, for some bounded positive integers $L, S$.
\item [$(ii)$] $r(\mat{X}^{l}_{\theta}) = \text{card}(\mathcal{I}_l)$ and $r(\mat{Y}^{s}_{\theta}) = \text{card}(\mathcal{J}_s)$.
\item [$(iii)$] $\norm{\mat{X}^{l}_{\theta}} < e^{-l + 1}$ and $\norm{\mat{Y}^{s}_{\theta}} < e^{-s + 1}$ by definition of sets $\mathcal{I}_l$ and $\mathcal{J}_s$ respectively in \eqref{eq:def_sets1}-\eqref{eq:def_sets2}.
\end{itemize}
In our analysis, we select $L$ and $S$ such that $\norm{\mat{X}_{\theta}^c} \leq \varepsilon$ and $\norm{\mat{Y}_{\theta}^c} \leq \varepsilon$. By definition of $\mathcal{I}_{l}, ~\mathcal{J}_s$, one can deduce that $\norm{\mat{X}_{\theta}^c} \leq e^{L+1}$ and $\norm{\mat{Y}_{\theta}^c} \leq e^{S+1}$, for fixed $L$ and $S$. To this end, we require $L = \left \lceil \ln\left(\frac{e}{\varepsilon}\right) \right \rceil$ and $S = \left \lceil \ln\left(\frac{e}{\varepsilon}\right) \right \rceil$ for both $\norm{\mat{X}_{\theta}^c} \leq \varepsilon$ and $\norm{\mat{Y}_{\theta}^c} \leq \varepsilon$ to be satisfied.

Based on the definitions above and using triangle inequality on $\norm{\mat{X}_{\theta} \left(\mat{G}^\top \mat{G} - \mat{I}\right)\mat{Y}_{\theta}} $, we obtain:
\begin{align}
\norm{\mat{X}_{\theta} \left(\mat{G}^\top \mat{G} - \mat{I}\right)\mat{Y}_{\theta}} &\stackrel{(i)}{=} \norm{\left(\sum_{l = 1}^{L} \mat{X}^{l}_{\theta}\right) \left(\mat{G}^\top \mat{G} - \mat{I}\right)\mat{Y}_{\theta} } \nonumber \\
&\stackrel{(iv)}{\leq} \sum_{l = 1}^{L} \norm{\mat{X}^{l}_{\theta}\left(\mat{G}^\top \mat{G} - \mat{I}\right)\mat{Y}_{\theta}} \nonumber \\ 
&\stackrel{(i)}{=} \sum_{l = 1}^{L}\norm{\mat{X}^{l}_{\theta} \left(\mat{G}^\top \mat{G} - \mat{I}\right)\left(\sum_{s = 1}^{S} \mat{Y}^{s}_{\theta}\right)} \nonumber \\
&\stackrel{(iv)}{\leq} \sum_{s = 1}^{S} \left(\sum_{l = 1}^{L} \norm{\mat{X}^{l}_{\theta}\left(\mat{G}^\top \mat{G} - \mat{I}\right)\mat{Y}^{s}_{\theta}}\right). \label{eq:result_2}
\end{align} 

For a pair $(l,~s)$ of positive integers, we define the following event, over the probability space defined by $\mat{G}$:
\begin{small}
\begin{align}
&\Omega(l,~s) \nonumber \\ &:= \left\{ \norm{\mat{X}^{l}_{\theta}\left(\mat{G}^\top \mat{G} - \mat{I}\right) \mat{Y}^{s}_{\theta}} \geq \left( \varepsilon + c_2\frac{r(\mat{X}^{l}_{\theta}) + r(\mat{Y}^{s}_{\theta})}{ c_1 \varepsilon t} \right) \norm{\mat{X}^{l}_{\theta}} \norm{\mat{Y}^{s}_{\theta}}\right\}, \nonumber
\end{align}
\end{small} where $c_1,~c_2 > 0$ are positive constants.
According to Lemma~\ref{lem:mm:rank}, the above event holds with probability:
\begin{align}
\mathbb{P}\left(\Omega(l, s)\right) &\leq c_2^{r(\mat{X}_{\theta}^{l}) + r(\mat{Y}_{\theta}^s)}\cdot e^{-c_1 \gamma_{l,s}t \cdot \min(\gamma_{l,s}, 1)} \label{eq:proba1}
\end{align}
where $\gamma_{l,s} = \varepsilon + c_2\frac{r(\mat{X}^{l}_{\theta}) + r(\mat{Y}^{s}_{\theta})}{ c_1 \varepsilon t}$. A useful observation for \eqref{eq:proba1} is given in the next lemma; the proof is provided in the appendix.
\begin{lemma}{\label{lem:prob_bound}}
Fix integer $t > 0$. The following inequality holds
\begin{align}{\label{eq:proba3}}
\Prob{\Omega(l, s)} &\leq e^{-c_1 t \cdot \min(\varepsilon, \varepsilon^2)}, \quad \forall l, s.
\end{align}
\end{lemma}

\begin{figure*}[tb]
\begin{align}
\|\mat{X}_{\theta} \big(\mat{G}^\top \mat{G} &- \mat{I}\big)\mat{Y}_{\theta}\|_2 \leq  \sum_{s = 1}^{S} \sum_{l = 1}^{L} \left( \left( \varepsilon + c_2\frac{r(\mat{X}^{l}_{\theta}) + r(\mat{Y}^{s}_{\theta})}{ c_1 t} \right) \norm{\mat{X}^{l}_{\theta}} \norm{\mat{Y}^{s}_{\theta}} \right) \nonumber \\ 
&= \varepsilon \sum_{s = 1}^{S} \sum_{l = 1}^{L} \norm{\mat{X}^{l}_{\theta}} \norm{\mat{Y}^{s}_{\theta}}  + \frac{c_2}{c_1 t}\sum_{s = 1}^{S} \sum_{l = 1}^{L} r(\mat{X}^{l}_{\theta})\norm{\mat{X}^{l}_{\theta}} \norm{\mat{Y}^{s}_{\theta}} + \frac{c_2}{c_1 t}\sum_{s = 1}^{S} \sum_{l = 1}^{L} r(\mat{Y}^{s}_{\theta}) \norm{\mat{X}^{l}_{\theta}} \norm{\mat{Y}^{s}_{\theta}} \nonumber \\
&\stackrel{(iii)}{<} \varepsilon\left(\sum_{s = 1}^{S} e^{-s + 1}\right)\left(\sum_{l = 1}^{L} e^{- l + 1}\right) + \frac{c_2}{c_1 t} \left(\sum_{s = 1}^{S} e^{-s + 1}\right) \left(\sum_{l = 1}^{L} r(\mat{X}^{l}_{\theta})\norm{\mat{X}^{l}_{\theta}}\right) + \frac{c_2}{c_1 t} \left(\sum_{l = 1}^{L} e^{-l + 1}\right) \left(\sum_{s = 1}^{S}r(\mat{Y}^{s}_{\theta}) \norm{\mat{Y}^{s}_{\theta}}\right) \nonumber \\ 
&\leq \varepsilon\left(\sum_{s = 1}^{\infty} e^{-s + 1}\right)\left(\sum_{l = 1}^{\infty} e^{- l + 1}\right) + \frac{c_2}{c_1 t} \left(\sum_{s = 1}^{\infty} e^{-s + 1}\right) \left(\sum_{l = 1}^{L} r(\mat{X}^{l}_{\theta})\norm{\mat{X}^{l}_{\theta}}\right) + \frac{c_2}{c_1 t} \left(\sum_{l = 1}^{\infty} e^{-l + 1}\right) \left(\sum_{s = 1}^{S}r(\mat{Y}^{s}_{\theta}) \norm{\mat{Y}^{s}_{\theta}}\right) \nonumber \\ 
&= \frac{\varepsilon \cdot e^2}{(e - 1)^2} + \frac{c_2}{c_1 t}\cdot\frac{e}{e - 1} \left[\left(\sum_{l = 1}^{L} r(\mat{X}^{l}_{\theta})\norm{\mat{X}^{l}_{\theta}}\right) + \left(\sum_{s = 1}^{S}r(\mat{Y}^{s}_{\theta}) \norm{\mat{Y}^{s}_{\theta}}\right)\right] \label{eq:result_3}
\end{align} 
\hrulefill
\end{figure*}
By conditioning on the event $ \bigcap_{s = 1}^{S}\bigcap_{l = 1}^{L} \Omega(l,s)^c $, we can upper bound \eqref{eq:result_2} as in \eqref{eq:result_3} where the infinite series $\sum_{i = 1}^{\infty} e^{ - i+1} = \frac{e}{e- 1}$ is used in the last equality.

To proceed, we observe the following for $\sum_{l = 1}^{L} r(\mat{X}^{l}_{\theta})\norm{\mat{X}^{l}_{\theta}}$; similar reasoning applies for $\sum_{s = 1}^{S}r(\mat{Y}^{s}_{\theta}) \norm{\mat{Y}^{s}_{\theta}}$. By construction in $(iii)$: 
\begin{align}
\norm{\mat{X}^{l}_{\theta}} < e^{-l+1} = e \cdot e^{-l} \leq e \cdot \min_{i\in \mathcal{I}_l} \sigma_i. \nonumber
\end{align} Thus, it is obvious that:
\begin{align}
r(\mat{X}^{l}_{\theta}) \cdot \norm{\mat{X}^{l}_{\theta}} < e \cdot r(\mat{X}^{l}_{\theta}) \cdot \min_{i\in \mathcal{I}_l} \sigma_i \leq e\cdot \|\mat{X}^{l}_{\theta}\|_{\star}, \nonumber
\end{align} since $\|\mat{E}\|_{\star} \geq  r(\mat{E}) \cdot \min_j \sigma_j(\mat{E})$, for any matrix $\mat{E}$. To this end,
\begin{align}
\sum_{l = 1}^{L} r(\mat{X}^{l}_{\theta})\norm{\mat{X}^{l}_{\theta}} < e \sum_{l = 1}^{L} \|\mat{X}^{l}_{\theta}\|_{\star} = e \|\mat{X}_{\theta}\|_{\star}.  \nonumber
\end{align} Similarly, we have $\sum_{s = 1}^{S}r(\mat{Y}^{s}_{\theta}) \norm{\mat{Y}^{s}_{\theta}} < e \|\mat{Y}_{\theta}\|_{\star}$ and, therefore, \eqref{eq:result_3} becomes:
\begin{align}
\norm{\mat{X}_{\theta} \left(\mat{G}^\top \mat{G} - \mat{I}\right)\mat{Y}_{\theta}} &< \frac{\varepsilon \cdot e^2}{(e - 1)^2} + \frac{c_2}{c_1 t}\cdot\frac{e}{e - 1} \left(\|\mat{X}\|_{\star} + \|\mat{Y}\|_{\star}\right) \nonumber \\ 
&= \frac{\varepsilon \cdot e^2}{(e - 1)^2} + \frac{c_2}{c_1 t}\cdot\frac{e}{e - 1} \left(\nr{\mat{X}} + \nr{\mat{Y}}\right) \nonumber \\ &\leq 3 \varepsilon \label{eq:result_4}
\end{align} where the last equality is satisfied since $t\geq \theta$ ($c_3\geq \frac{c_2e}{c_1 (e - 1)}$). By the union bound and since $0 < \varepsilon < 1$, \eqref{eq:result_4} is violated with probability:
\begin{align}
\mathbb{P}\left(\left(\bigcap_{s = 1}^{S}\bigcap_{l = 1}^{L} \Omega(l,s)^c \right)^c\right) &\stackrel{\eqref{eq:proba3}}{\leq} \sum_{s = 1}^{S} \sum_{l = 1}^{L} e^{-c_1 \cdot \varepsilon^2 \cdot t} \nonumber \\ &= e^{-c_1 \cdot \varepsilon^2 \cdot t + \ln\left( \left\lceil \ln\left(\frac{e}{\varepsilon}\right)  \right\rceil \right) } \leq \frac{\delta}{2} \nonumber
\end{align} where the last inequality is satisfied since $t \geq \frac{\ln(2/\delta) + \ln\left( \left\lceil \ln(e/\varepsilon)\right\rceil\right)}{c_1 \cdot \varepsilon^2}$.

The terms in \eqref{eq:termConstant} can be bounded as follows:
\begin{align}
\norm{\mat{X}_{\theta}^c \mat{Y}_{\theta}} &+ \norm{\mat{X}_{\theta} \mat{Y}_{\theta}^c} + \norm{\mat{X}_{\theta}^c \mat{Y}_{\theta}^c} \nonumber \\ 
&\leq \norm{\mat{X}_{\theta}^c} \norm{\mat{Y}_{\theta}} + \norm{\mat{X}_{\theta}} \norm{\mat{Y}_{\theta}^c} + \norm{\mat{X}_{\theta}^c} \norm{\mat{Y}_{\theta}^c} \nonumber \\
&\leq 2\varepsilon + \varepsilon^2 \leq 3 \varepsilon \nonumber
\end{align}
where we used the fact that $\norm{\mat{X}_\theta},~\norm{\mat{Y}_\theta} \leq 1$ and the inequalities $\norm{\mat{X}_\theta^c} \leq \varepsilon$ (similarly for $\mat{Y}_\theta^c$).

To bound the terms in \eqref{eq:termNormGaussian}, we provide the next Lemma; the proof is given in the appendix.
\begin{lemma}\label{lem:tailTerms}
If $t\geq \theta + 8 \ln(8/\delta)$ and $0< \varepsilon <1 $, then, with probability at least $\delta / 2$, the term appearing in \eqref{eq:termNormGaussian} is at most $33\varepsilon$.
\end{lemma}
Applying the union bound on Lemma~\ref{lem:tailTerms} and the complement of $ \bigcap_{s = 1}^{S}\bigcap_{l = 1}^{L} \Omega(l,s)^c $, we conclude that the following inequality holds:
\begin{align}
&\norm{\mat{X} \mat{G}^\top \mat{G} \mat{Y} - \mat{X} \mat{Y}} \leq 39 \varepsilon, \nonumber
\end{align} with probability at least $1 - \delta$. By rescaling $\varepsilon$ we obtain the required result.
