mathx”17

# Robust recovery of low-rank matrices and low-tubal-rank tensors from noisy
sketches††thanks: A.M. was supported by U.S. Air Force Grant FA9550-18-1-0031,
NSF Grant DMS 2027299, and the NSF+Simons Research Collaborations on the
Mathematical and Scientific Foundations of Deep Learning. Y.Z. was partially
supported by NSF-Simons Research Collaborations on the Mathematical and
Scientific Foundations of Deep Learning. Y.Z. also acknowledges support from
NSF DMS-1928930 during his participation in the program “Universality and
Integrability in Random Matrix Theory and Interacting Particle Systems” hosted
by the Mathematical Sciences Research Institute in Berkeley, California,
during the Fall semester of 2021.

Anna Ma Department of Mathematics, University of California Irvine
(anna.ma@uci.edu). Dominik Stöger Department of Mathematics and Mathematical
Institute for Machine Learning and Data Science (MIDS), KU Eichstätt-
Ingolstadt (dominik.stoeger@ku.de). Yizhe Zhu Department of Mathematics,
University of California Irvine (yizhe.zhu@uci.edu).

###### Abstract

A common approach for compressing large-scale data is through matrix
sketching. In this work, we consider the problem of recovering low-rank
matrices from two noisy linear sketches using the double sketching algorithm
discussed in Fazel et al. (2008). Using tools from non-asymptotic random
matrix theory, we provide the first theoretical guarantees characterizing the
error between the output of the double sketch algorithm and the ground truth
low-rank matrix. We apply our result to the problems of low-rank matrix
approximation and low-tubal-rank tensor recovery.

##  1 Introduction

The prevalence of large-scale data in data science applications has created
immense demand for methods that can aid in reducing the computational cost of
processing and storing said data. Oftentimes, data, such as images, videos,
and text documents, can be represented as a matrix, and thus, the ability to
store matrices efficiently becomes an important task. One way to efficiently
store a large-scale matrix $X_{0}$ $\in\mathbb{R}^{m\times n}$ is to store a
sketch of the matrix, i.e., another matrix $Y\in\mathbb{R}^{m_{1}\times
n_{1}}$ such that two goals are accomplished. Firstly, the sketch of $X_{0}$
must be cheaper to store than $X_{0}$ itself, i.e., we want $m_{1}n_{1}\ll
mn$. Second, we must be able to recover $X_{0}$ from its sketch, meaning that
we are able to reconstruct $X_{0}$ exactly or at least be able to obtain a
high-quality approximation of $X_{0}$.

A variety of works have been produced for the setting in which $X_{0}$ is a
low-rank matrix, and one wishes to recover $X_{0}$ using its sketches, see
e.g. [16, 46, 11, 34, 38, 39]. However, in certain settings, one may only have
access to noisy sketches. For example, noise in the sketch can arise from
errors in the sensors or quantization, a low-rank approximation of the ground
truth matrix, or be a result of sketching noisy matrices [45, 6]. Noisy
sketches have also appeared in recovering low-rank matrices from noisy linear
observations [29, 35]. It should be noted that in the aforementioned noisy
observation models, only one noisy sketch is obtained, while in the model we
consider here (4), there are two noisy sketches available. Noisy sketches are
also used in differentially private algorithms for matrix factorization [40,
1]. In particular, the model we consider in (4) is closely related to the
output perturbation method described in [40, 1], which is an important
intermediate step for algorithms that preserve differential privacy of a low-
rank matrix factorization (see Algorithm 1 in [40]).

In this work, we analyze the noisy double-sketch algorithm originally proposed
but not theoretically studied in [16]. Here, the term “double-sketch” refers
to the first algorithm in [16] where two independent random sketches matrices
are applied to $X_{0}$ by multiplying from the left and right, respectively.
We show that when the sketching matrices are i.i.d. complex Gaussian random
matrices, one can recover the original low-rank matrix $X_{0}$ with high
probability where the error on the approximation depends on the noise level
for both sketches. Here, we do not assume that one has access to the exact
rank of $X_{0}$ but instead an approximate rank $r$, which satisfies
$r>r_{0}:=\text{rank}(X_{0})$. We also remark on the utility of our
theoretical guarantees for when the double sketch algorithm is used not for
low-rank matrix recovery but instead for low-rank approximation with noise.
Lastly, we present results for the application of this work to a more extreme
large-scale data setting in which one wants to recover a low-tubal-rank tensor
[23] (see Definition 5), which are used as a model for multidimensional data
in applications.

A key step in our robust recovery analysis is to control the perturbation
error of a low-rank matrix under noise. A standard way is to apply Wedin’s
Theorem, or Davis-Kahan Theorem [44, 12, 10, 30], which results in a bound
that depends on the condition number of the low-rank matrix. However, our
proof is based on an exact formula to calculate the difference between output
and the ground truth matrix and a detailed analysis of the random matrices
(extreme singular value bounds of Gaussian matrices [42, 33, 37] and the least
singular value of truncated Haar unitary matrices [5, 13]) involved in the
double sketch algorithm. This novel approach yields a bound independent of the
condition number of the low-rank matrix (Theorem 2). Due to the Gaussian
structure of our sensing matrices, our results are non-asymptotic, and all the
constants involved in the probabilistic error bounds are explicit.

###  1.1 Low-rank matrix reconstruction

The double sketch method, which is analyzed in this paper, was proposed (but
not theoretically studied) in [16] to recover low-rank matrices from noisy
sketches. This method is a single-pass method, meaning that the input matrix
needs only to be read once. Related single-pass sketching methods were
analyzed in [11, 46, 38, 39] in the noise-free scenario, meaning that, in
contrast to the model studied in this paper, the sketches are not perturbed by
additive noise.

The single-pass model is motivated by the fact that in many modern
applications, one cannot make several passes over the input matrix $X_{0}$.
However, if one is able to make several passes over the input matrix $X_{0}$,
one can, in many applications, expect superior performance compared to single-
pass algorithms. We note that multi-pass algorithms for low-rank matrix
approximation using sketches have been studied in, e.g., [51, 19, 18].

The so-called problem of compressive PCA was studied in [34] and [2]. It can
be interpreted as a variant of sketching where only the columns of a matrix
are sketched. However, this problem is not directly comparable to the setting
in the paper at hand, as in compressive PCA, a different sketching matrix is
used for each column.

###  1.2 Low-tubal-rank tensor recovery

The notion of a low-tubal-rank tensor stems from the t-product, originally
introduced by [24]. We state the relevant definitions for order-3 tensors. We
refer to [23] for a more comprehensive treatment. More general definitions for
tensors of higher orders can be found in [22, 25].

######  Definition 1 (Operations on tensors, see, e.g., [24, Section 3.2]).

Let $\mathcal{A}\in\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$. The unfold of
a tensor is defined to be the frontal slice stacking of that tensor. In other
words,

$\text{unfold}({\mathcal{A}})=\begin{pmatrix}A_{1}\\\ A_{2}\\\ \vdots\\\
A_{n_{3}}\\\ \end{pmatrix}\in\mathbb{C}^{n_{1}n_{3}\times n_{2}},$

where $A_{i}=\mathcal{A}_{:,:,i}$ denotes the $i^{th}$ frontal slice of
$\mathcal{A}$. We define the inverse of the $\text{unfold}(\cdot)$ as
$\text{fold}{(\cdot)}$ so that
$\text{fold}(\text{unfold}(\mathcal{A}))=\mathcal{A}$. The block circulant
matrix of $\mathcal{A}$ is:

$\text{bcirc}({\mathcal{A}})=\begin{pmatrix}A_{1}&A_{n_{3}}&A_{n_{3}-1}&\ldots&A_{2}\\\
A_{2}&A_{1}&A_{n_{3}}&\ldots&A_{3}\\\ \vdots&\vdots&\vdots&\ddots&\vdots\\\
A_{n_{3}}&A_{n_{3}-1}&A_{n_{3}-2}&\ldots&A_{1}\\\
\end{pmatrix}\in\mathbb{C}^{n_{1}n_{3}\times n_{2}n_{3}}.$

The conjugate transpose of a tensor $\mathcal{A}\in\mathbb{C}^{n_{1}\times
n_{2}\times n_{3}}$ is the tensor $\mathcal{A}^{*}\in\mathbb{C}^{n_{2}\times
n_{1}\times n_{3}}$ obtained by conjugate transposing each of the frontal
slices and then reversing the order of transposed frontal slices $2$ through
$n_{3}$.

######  Definition 2 (Tensor t-product, see [24, Section 3.2]).

Let $\mathcal{A}\in\mathbb{C}^{n_{1}\times\ell\times n_{3}}$ and
$\mathcal{B}\in\mathbb{C}^{\ell\times n_{2}\times n_{3}}$ then the t-product
between $\mathcal{A}$ and $\mathcal{B}$, denoted $\mathcal{A}*\mathcal{B}$, is
a tensor of size $n_{1}\times n_{2}\times n_{3}$ and is computed as:

$\mathcal{A}*\mathcal{B}=\text{fold}(\text{bcirc}(\mathcal{A})\text{unfold}(\mathcal{B})).$
(1)

######  Definition 3 (Mode-3 fast Fourier transformation (FFT), see [23,
Section 2]).

The mode-3 FFT of a tensor $\mathcal{A}$, denoted $\widehat{\mathcal{A}}$, is
obtained by applying the discrete Fourier Transform matrix,
$F\in\mathbb{C}^{n_{3}\times n_{3}}$, to each $\mathcal{A}_{i,j,:}$ of
$\mathcal{A}$:

$\widehat{\mathcal{A}}_{i,j,:}=F\mathcal{A}_{i,j,:}.$ (2)

Here, $F$ is a unitary matrix, $\mathcal{A}_{i,j,:}$ is an $n_{3}$-dimensional
vector, and the product is the usual matrix-vector product. Similarly, the
inverse mode-3 fast Fourier transform is defined by applying the inverse
discrete Fourier Transform matrix to each $\widehat{\mathcal{A}}_{i,j,:}$ of
$\widehat{\mathcal{A}}$.

######  Definition 4 (t-SVD, see [24, Theorem 4.1]).

The Tensor Singular Value Decomposition (t-SVD) of a tensor
$\mathcal{M}\in\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$ is given by

$\mathcal{M}=\mathcal{U}*\mathcal{S}*\mathcal{V}^{*},$ (3)

where $\mathcal{U}\in\mathbb{C}^{n_{1}\times n_{1}\times n_{3}}$ and
$\mathcal{V}\in\mathbb{C}^{n_{2}\times n_{2}\times n_{3}}$ are unitary tensors
and $\mathcal{S}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ is a tubal
tensor (a tensor in which each frontal slice is diagonal), and $*$ denotes the
t-product.

######  Definition 5 (Tubal rank, see [23, Section 4.1]).

The tubal rank of a tensor
$\mathcal{M}=\mathcal{U}*\mathcal{S}*\mathcal{V}^{*}$ is the number of non-
zero singular tubes of $\mathcal{S}$, and it is uniquely defined.

######  Definition 6 (CP rank).

The CP rank of an order three tensor $\mathcal{M}$ is the smallest integer $r$
such that $\mathcal{M}$ is a sum of $r$ rank-1 tensors:

$\mathcal{M}=\sum_{i=1}^{r}u_{i}\otimes v_{i}\otimes w_{i},$

where
$u_{i}\in\mathbb{C}^{n_{1}},v_{i}\in\mathbb{C}^{n_{2}},w_{i}\in\mathbb{C}^{n_{3}}$,
$1\leq i\leq r$.

If a tensor $\mathcal{M}$ has CP rank $r$ then its tubal rank is at most $r$,
see [50, Remark 2.3].

######  Definition 7 (Tensor Frobenius norm, see, e.g., [24, Definition 2.1]).

Let $\mathcal{M}\in\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$. The Frobenius
norm of $\mathcal{M}$ is given by

$\|\mathcal{M}\|_{F}=\left(\sum_{i_{1},i_{2},i_{3}}|\mathcal{M}_{i_{1},i_{2},i_{3}}|^{2}\right)^{1/2}.$

Other low-rank tensor sketching approaches have been proposed for low-CP-rank
tensors [21] and low-Tucker-rank tensors [36]. In the following, we focus on
low-tubal-rank tensors since this is the topic of this paper.

Related work considered the recovery of low-tubal-rank tensors through general
linear Gaussian measurements of the form $y=A\text{vec}(\mathcal{X})$ [26, 49,
48]. This can be seen as a generalization of the low-rank matrix recovery
problem [32] to low-tubal-rank tensors. The proof of tensor recovery under
Gaussian measurements in [26, 49, 48] relies crucially on the assumption that
the entries of the measurement matrix $A$ are i.i.d. Gaussian. In this
setting, it was shown that the tensor nuclear norm is an atomic norm, and a
general theorem from [9, Corollary 12] for i.i.d. measurements for atomic
norms was used to establish recovery guarantees. In [43], a non-convex
surrogate for the tensor nuclear norm was proposed and studied.

An extension of the matrix sketching algorithm in [39] to a low-tubal-rank
approximation of tensors was considered in [31]. However, their setting does
not cover noisy sketching, which is the topic of this paper. Streaming low-
tubal-rank tensor approximation was considered in [47].

#### Notations

We define a standard complex Gaussian random variable $Z$ as $Z=X+iY$, where
$X\sim N(0,1/2),Y\sim N(0,1/2)$ and $X,Y$ are independent. By $\sigma_{i}(A)$,
we define the $i$-th largest singular value of a matrix $A$. By
$\sigma_{\min}(A)$, we denote the smallest non-zero singular value of a matrix
$A$. $\|A\|_{2\rightarrow 2}$ is the spectral norm of a matrix $A$, and
$\vvvert A\vvvert$ is a general norm of $A$. The matrix $A^{*}$ is the complex
conjugate, and $A^{\dagger}$ is the pseudo-inverse of $A$. Let $U$ be a matrix
with orthonormal columns. By $U_{\perp}$, we denote the orthogonal complement
of $U$, which means the column vectors of $U$ and the column vectors of
$U_{\perp}$ form a complete orthonormal basis.

#### Organization of the paper

The rest of the paper is organized as follows. We provide more background on
low-rank matrix recovery from sketches in Section 2. We state our main results
and corollaries of robust recovery of low-rank matrix and low-tubal-rank
tensors in Section 3. The proofs of all theorems and corollaries are provided
in Section 4, and auxiliary lemmas are provided in the Appendix.

##  2 Problem formulation and background

In the following, we denote by $X_{0}\in\mathbb{C}^{n_{1}\times n_{2}}$ the
target matrix with rank $r_{0}$ to be sketched. For that, we use two sketching
matrices $S\in\mathbb{C}^{r\times n_{1}}$ and $\tilde{S}\in\mathbb{C}^{r\times
n_{2}}$, which are assumed to be two independent complex Gaussian random
matrices with $r\geq r_{0}$. We consider the model

$\displaystyle Y:=SX_{0}+Z,\quad\tilde{Y}:=\tilde{S}X_{0}^{*}+\tilde{Z},$ (4)

where $Z\in\mathbb{C}^{r\times n_{2}}$ and $\tilde{Z}\in\mathbb{C}^{r\times
n_{1}}$ represent additive noise. This model has first been considered (but
not been analyzed) in [16]. The matrices $SX_{0}$ and $\tilde{S}X_{0}^{*}$,
represent the noiseless sketches, which contain information about the row
space and column space of the target matrix $X_{0}$. Our goal is to estimate
the true matrix $X_{0}$ given the two noisy sketches $Y$ and $\tilde{Y}$. For
that purpose, we use the double-sketch method proposed in [16], which outputs

$X=\tilde{Y}^{*}(S\tilde{Y}^{*})^{\dagger}Y.$ (5)

The matrix $X$ represents an estimate of the target matrix $X_{0}$.

Note that formula (5) has also been studied in [11, Theorem 4.7] in the
context of low-rank approximation, i.e., $Z=0$, $\tilde{Z}=0$ and $X_{0}$ is
approximately low-rank. (See also Section 3.2, where we apply our results to
the setting of noisy low-rank matrix approximation.) In [38, Section 4], an
algebraically equivalent, but more numerical stable reformulation of equation
(5) was described, again for the problem of low-rank matrix approximation. We
outline the main ideas of this approach below. Moreover, a similar approach as
formula (5) has been studied in [46]. Namely, in that paper instead of (5) the
formula $X=\tilde{Y}_{r_{0}}^{*}(S\tilde{Y}_{r_{0}}^{*})^{\dagger}Y_{r_{0}}$
has been used, where by $Y_{r_{0}}$, respectively $\tilde{Y}_{r_{0}}$, we
denote the truncated $r_{0}$-term singular value decomposition of $Y$,
respectively $\tilde{Y}$.

In the noiseless case when $Z=0$, $\tilde{Z}=0$, we denote the output of (5)
as $\overline{X}$. In this case, the output will read as

$\displaystyle\overline{X}=X_{0}\tilde{S}^{*}(SX_{0}\tilde{S}^{*})^{\dagger}SX_{0}.$
(6)

While this work is concerned with the theoretical analysis of the double-
sketch method (5), we want to make a few comments regarding its practical
implementation. Our presentation will be based on [38], which we refer to for
more details. First of all, note that in many applications, the matrix $X_{0}$
is too large to be stored on a computer. A typical scenario, in which the
sketching method is used, is a streaming model (see [28]). That is, we have

$X_{0}=H_{1}+H_{2}+H_{3}+H_{4}+\ldots,$

where the updates $H_{i}$ must be discarded after they have been processed.
Typically the matrices $H_{i}$ possess a simple structure, such as being rank-
one. This means that in this scenario, we can obtain the sketches $SX_{0}$ and
$\tilde{S}X_{0}^{*}$ without ever needing to compute or store the full matrix
$X_{0}$. Note that the sketching method presented in this paper is
particularly suited for scenarios where we can only make one single pass over
the data as in the streaming model.

In order to compute $X$ from formula (5) one might obtain the matrix
$(S\tilde{Y}^{*})^{\dagger}Y$ via its associated least-squares problem.
However, the matrix $S\tilde{Y}^{*}$ might be poorly conditioned, leading to
an inaccurate solution. As suggested in [38], one can avoid this by
considering the unitary-triangular decomposition $\tilde{Y}^{*}=QR$, where
$Q\in\mathbb{C}^{n_{1}\times r}$. Then, (if $\tilde{Y}$ has full rank and
$SQ\in\mathbb{C}^{r\times r}$ is invertible) we observe that (5) is equivalent
to

$X=Q(SQ)^{\dagger}Y.$

Note that $SQ$ is typically better conditioned than $S\tilde{Y}^{*}$, implying
that $(SQ)^{\dagger}Y$ can be computed more accurately. The pseudo-code for
this method is provided in Algorithm 1.

Algorithm 1 Matrix Recovery from Double-Sketch

Sketches $Y$, $\tilde{Y}$, and sketching matrix $S$ from (4)

$X\in\mathbb{R}^{n_{1}\times n_{2}}$ approximation of target matrix $X_{0}$

function doubleSketchMatrixRecovery($Y$, $\tilde{Y}$, $S$)

Compute the QR Decomposition of $\tilde{Y}$: $\tilde{Y}^{*}=QR$ where
$Q\in\mathbb{R}^{n_{1}\times r}$ and $R\in\mathbb{R}^{r\times r}$

$X=Q(SQ)^{\dagger}Y$

return $X$

end function

In the end, let us stress that in order to access individual entries of $X$ it
suffices to compute the corresponding inner products between $Q$ and
$(SQ)^{\dagger}Y$. In particular, there is no need to compute or store the
complete matrix $X$.

##  3 Main results

###  3.1 Low-rank matrix recovery from noisy sketches

In our first result, we show that without noise (5) exactly recovers $X_{0}$,
i.e., $\overline{X}=X_{0}$ with probability 1.

######  Theorem 1 (Exact recovery).

Let $S\in\mathbb{C}^{r\times n_{1}},\tilde{S}\in\mathbb{C}^{r\times n_{2}}$ be
two independent complex standard Gaussian random matrices. Furthermore, let
$X_{0}\in\mathbb{C}^{n_{1}\times n_{2}}$ be a matrix with rank $r_{0}$. If
$r\geq r_{0}$ and $Z=0$, $\tilde{Z}=0$ then with probability one
$\overline{X}=X_{0}$, where $\overline{X}$ is as defined in (6).

Our Theorem 1 generalized the exact recovery result [16, Lemma 6], where $r$
is assumed to be exactly $r_{0}$. Our Theorem 1 implies the exact value of
$r_{0}$ is not needed for the double sketch algorithm, and one can always use
the parameter $r\geq r_{0}$. In fact, our robust recovery result (Theorem 2)
suggests choosing a larger $r$ makes the output of the double sketch algorithm
more robust to noise.

When $Z,\tilde{Z}$ are not all zero, the robust recovery guarantee is given as
follows.

######  Theorem 2 (Robust recovery).

Let $S\in\mathbb{C}^{r\times n_{1}}$, $\tilde{S}\in\mathbb{C}^{r\times n_{2}}$
be two independent standard complex Gaussian matrices. Let
$Z\in\mathbb{C}^{r\times n_{1}}$ be any matrix, and
$\tilde{Z}\in\mathbb{C}^{r\times n_{2}}$ be a matrix such that
$(\tilde{S},\tilde{Z})$ is independent of $S$, and
$\tilde{Y}=\tilde{S}X^{*}_{0}+\tilde{Z}$ is almost surely of rank $r$ with
$r_{0}<r<n_{1}$. For any $\delta_{1},\delta_{2},\epsilon>0$ such that
$1>\delta_{2}>\exp\left(-\left(\sqrt{r}-\sqrt{r_{0}}\right)^{2}\right)$ and
$\epsilon<1$, with probability at least $1-\delta_{1}-\delta_{2}-\epsilon$
over the randomness of $\tilde{S}$ and $S$, the output $X$ from Algorithm 1
satisfies

$\displaystyle\vvvert X-X_{0}\vvvert$
$\displaystyle\leq\frac{\sqrt{r(n_{1}-r)}\vvvert\tilde{Z}\vvvert}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{0}}-\sqrt{\log(1/\delta_{2})})}+\frac{\sqrt{r}\vvvert
Z\vvvert}{\sqrt{\log(1/(1-\epsilon))}},$

where $\vvvert\cdot\vvvert$ is any matrix norm that satisfies for any two
matrices $A$ and $B$, $\vvvert AB\vvvert\leq\|A\|_{2\rightarrow 2}\vvvert
B\vvvert$, and $\vvvert A\vvvert=\vvvert A^{*}\vvvert$. In particular, it
holds for $\|\cdot\|_{2\rightarrow 2}$ and $\|\cdot\|_{F}$.

The condition that $(\tilde{S},\tilde{Z})$ is independent of $S$ is important
in our proof to apply the least singular value bound for truncated Haar random
matrices, see Lemma 4. The condition that $\tilde{Y}$ is of rank $r$ can be
easily verified in different settings. For example, it holds when $\tilde{Z}$
is Gaussian noise independent of $\tilde{S}$, or
$\tilde{Y}=\tilde{S}(X_{0}^{*}+\hat{Z})$, where $X_{0}^{*}+\hat{Z}$ is of full
rank.

Our proof of Theorem 2 works for $r=r_{0}$ or $r=n_{1}$, but the error bounds
are slightly different.

  * •

When $r=r_{0}<n_{1}$, (18) is replaced by the estimate that
$\sigma_{\min}(\tilde{S}V_{0})\geq\sqrt{\log(1/(1-\delta_{2}))}r^{-1/2}$ with
probability $1-\delta_{2}$. This implies with probability at least
$1-\delta_{1}-\delta_{2}-\epsilon$,

$\displaystyle\vvvert X-X_{0}\vvvert$
$\displaystyle\leq\frac{r\sqrt{n_{1}-r}\vvvert\tilde{Z}\vvvert}{\sqrt{\delta_{1}(\log(1/(1-\delta_{2}))}}+\frac{\sqrt{r}\vvvert
Z\vvvert}{\sqrt{\log(1/(1-\epsilon))}}.$

  * •

When $r=n_{1}$, $U_{A}$ defined in (11) is a unitary matrix and $S$ is
invertible. We find $\tilde{X}=X_{0}$ in (14) and following the rest of the
proof, we obtain with probability at least $1-\epsilon$,

$\vvvert X-X_{0}\vvvert\leq\frac{\sqrt{n_{1}}\vvvert
Z\vvvert}{\sqrt{\log(1/(1-\epsilon))}}.$

Note that in the case when $r=n_{1}$, the error is independent of $\tilde{Z}$.
This is because in (5), when $S$ is invertible, from part (1) in Lemma 2, the
expression becomes

$X=\tilde{Y}^{*}{(\tilde{Y}^{*})}^{\dagger}S^{-1}Y=S^{-1}Y,$

which is independent of $\tilde{Z}$.

###  3.2 Low-rank matrix approximation from noisy sketches

Note that Theorem 2 uses the assumption that the rank of the ground-truth
matrix $X_{0}$ is exactly $r_{0}$. However, in most practical scenarios the
matrix $X_{0}$ is only approximately of low-rank, i.e., the remaining singular
values decay quickly to zero. Indeed, there is a large literature focusing on
single-pass sketching algorithms for low-rank approximation (see, e.g., [46,
38, 39]). However, to the best of our knowledge, none of these works consider
the scenario that there is additive noise on the sketched matrices.

In the following, we describe how one can utilize Theorem 2 to also obtain a
bound for the approximation error for the scenario that $X_{0}$ is
approximately low-rank and that there is additive noise. In fact, when $X_{0}$
is not low-rank, we can write $X_{0}=X_{1}+E$, where $X_{1}$ is the best
rank-$r_{1}$ approximation of $X_{0}$. Letting $r_{0}>r>r_{1}$, we can use the
noisy double sketch model in (4) to consider the sketches

$\displaystyle Y$ $\displaystyle=SX_{0}+Z=SX_{1}+(SE+Z),$
$\displaystyle\tilde{Y}$
$\displaystyle=\tilde{S}X_{0}^{*}+\tilde{Z}=\tilde{S}X_{1}^{*}+(\tilde{S}E^{*}+\tilde{Z}).$

By interpreting $SE+Z$ and $\tilde{S}E^{*}+\tilde{Z}$ as perturbations of the
“true” sketches $SX_{1}$ and $\tilde{S}X_{1}^{*}$, we can apply Theorem 2,
which yields Corollary 1. The proof of Corollary 1 is given in Section 4.3.

######  Corollary 1 (Low-rank approximation with noisy sketches).

For $X_{0}\in\mathbb{C}^{n_{1}\times n_{2}}$ and an integer $r_{1}$ with that
$r_{1}<r<n_{1}$, consider the sketches

$Y=SX_{0}+Z,\quad\tilde{Y}=\tilde{S}X_{0}^{*}+\tilde{Z}.$

where $S$ and $\tilde{S}$ are two independent complex Gaussian random
matrices. Moreover, assume that $S$ is independent of $\tilde{Z}$. Recall that
the output of Algorithm 1 is given by

$X=\tilde{Y}^{*}(S\tilde{Y}^{*})^{\dagger}Y.$

For any $\delta_{1},\delta_{2},\epsilon>0$ such that
$1\geq\delta_{2}>\exp(-(\sqrt{r}-\sqrt{r_{1}})^{2})$ and $\epsilon<1$, with
probability at least $1-\delta_{1}-3\delta_{2}-\epsilon$ and when $\tilde{Y}$
is of rank $r$, the output $X$ satisfies

$\displaystyle\|X-X_{0}\|_{2\rightarrow
2}\leq\sigma_{r_{1}+1}\left(X_{0}\right)\cdot$
$\displaystyle\left(\frac{\sqrt{r(n_{1}-r)}(\sqrt{r}+\sqrt{n_{2}}+\sqrt{\log(1/\delta_{2})})}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})}+\frac{\sqrt{r}\left(\sqrt{r}+\sqrt{n_{1}}+\sqrt{\log(1/\delta_{2})}\right)}{\sqrt{\log(1/(1-\epsilon))}}+1\right)$
$\displaystyle+\frac{\sqrt{r(n_{1}-r)}\|\tilde{Z}\|_{2\rightarrow
2}}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})}+\frac{\sqrt{r}\|Z\|_{2\rightarrow
2}}{\sqrt{\log(1/(1-\epsilon))}}.$

######  Remark 1.

When $Z=0$ and $\tilde{Z}=0$ and $X_{0}$ is approximately low-rank, in [46]
and [38], two different algorithms have been proposed to approximate the low-
rank matrix $X_{0}$ using the two sketches $Y$ and $\tilde{Y}$. In the
noiseless setting, Corollary 1 below yields a weaker error bound when compared
to the bounds in [46] and [38]. However, our result can handle noise in the
two sketches $Y$ and $\tilde{Y}$, while the proofs from [46, 38] are not
applicable. In fact, the proofs in [46, 38] heavily rely on the assumption
$Z=0$ and $\tilde{Z}=0$ to use properties of orthogonal projections, which
only hold in this noiseless scenario. See for example [38, Fact A.2].

The error bound in Corollary 1 consists of two parts. The term involving
$\sigma_{r_{1}+1}\left(X_{0}\right)$ comes from comparing $X$ with the best
rank-$r_{1}$ approximation to $X_{0}$, and the term without
$\sigma_{r_{1}+1}\left(X_{0}\right)$ is the error from the noisy sketches.
When the rank of $X_{0}$ is $r_{1}$, the first error term is $0$, and
Corollary 1 reduces to Theorem 2. The bound in Corollary 1 is true for any
$r_{1}<r$. Therefore one can optimize $r_{1}$ to find the best bound in terms
of the failure probability and the approximation error.

###  3.3 Application to sketching low-tubal-rank tensors

The approach set forth in (5) can be used to sketch and recover low-tubal-rank
tensors. For such an application, one considers the low-tubal-rank tensor
$\mathcal{X}_{0}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ with tubal rank
$r_{0}$. Taking the mode-3 FFT of $\mathcal{X}_{0}$, one obtains
$\hat{\mathcal{X}}_{0}$ which is composed of a collection of $n_{3}$ matrices
(frontal slices) of dimension $n_{1}\times n_{2}$ with rank at most $r_{0}$.
As such, (5) can be used to sketch each of the $n_{3}$ frontal slices of
$\mathcal{X}_{0}$. The recovery of target tensor $\mathcal{X}_{0}$ from double
sketches (7) is outlined in Algorithm 2 and Corollary 2 presents guarantees
for the approximation error. The proof is given in Section 4.3.

Algorithm 2 Tensor Recovery from Double-Sketch

Sketches $\mathcal{Y}$, $\widetilde{\mathcal{Y}}$, and sketching tensor
$\mathcal{S}$ from (7)

$\mathcal{X}\in\mathbb{R}^{n_{1}\times n_{2}\times n_{3}}$ approximation of
target tensor $\mathcal{X}_{0}$

function doubleSketchTensorRecovery($\mathcal{Y}$, $\widetilde{\mathcal{Y}}$,
$\mathcal{S}$)

Compute $\hat{\mathcal{Y}}$, $\widehat{\widetilde{\mathcal{Y}}}$,
$\hat{\mathcal{S}}$ via mode-3 FFT $\triangleright$ See Definition 3.

for $i\in\\{1,\cdots,n_{3}\\}$ do

$\widehat{\mathcal{X}}_{i}=\textsc{DoubleSketchMatrixRecovery}(\widehat{\mathcal{Y}}_{i},\widehat{\widetilde{\mathcal{Y}}}_{i},\widehat{\mathcal{S}}_{i})$
$\triangleright$ See Algorithm 1.

end for

Compute $\mathcal{X}$ from $\widehat{\mathcal{X}}$ via the mode-3 inverse FFT

return $\mathcal{X}$

end function

######  Corollary 2 (Recovering low tubal-rank tensors).

Let $\mathcal{X}_{0}\in\mathbb{C}^{n_{1}\times n_{2}\times n_{3}}$ be a low-
tubal-rank tensor with rank $r_{0}$. Furthermore, let $S\in\mathbb{C}^{r\times
n_{1}},\tilde{S}\in\mathbb{C}^{r\times n_{2}}$ be two independent complex
standard Gaussian random matrices and $\mathcal{Z}\in\mathbb{R}^{r\times
n_{2}\times n_{3}}$ and $\mathcal{\tilde{Z}}\in\mathbb{R}^{r\times n_{1}\times
n_{3}}$ be two noise tensors such that $\mathcal{\tilde{Z}}$ is stochastically
independent of $S$. Consider the measurements

$\displaystyle\mathcal{Y}$
$\displaystyle=\mathcal{S}*\mathcal{X}_{0}+\mathcal{Z},\,\,\,\mathcal{\tilde{Y}}=\mathcal{\tilde{S}}*\mathcal{X}_{0}^{*}+\mathcal{\tilde{Z}},$
(7)

where $\mathcal{S}_{1}=S$, $\mathcal{\tilde{S}}_{1}=\tilde{S}$ and
$\mathcal{S}_{k}=\textbf{0}$, $\mathcal{\tilde{S}}_{k}=\textbf{0}$ for all
$k\in\\{2,\dots,n_{3}\\}$.

  1. 1.

(Exact recovery) If $\mathcal{Z}=\textbf{0}$ and
$\mathcal{\tilde{Z}}=\textbf{0}$ then with probability 1,
$\mathcal{X}=\mathcal{X}_{0}$.

  2. 2.

(Robust recovery) If $r_{0}<r<n_{1}$ and for all $k\in[n_{3}]$,
$\hat{\tilde{\mathcal{Y}}}_{k}$ is of rank $r$, then for any
$\delta_{1},\delta_{2},\epsilon>0$ such that
$1>\delta_{2}>\exp(-(\sqrt{r}-\sqrt{r_{0}})^{2})$ and $\epsilon<1$, with
probability at least $1-(\delta_{1}+\delta_{2}+\epsilon)n_{3}$,

$\displaystyle\|\mathcal{X}-\mathcal{X}_{0}\|_{F}^{2}\leq\frac{2{r(n_{1}-r)}\|{\tilde{\mathcal{Z}}}\|_{F}^{2}}{{\delta_{1}}(\sqrt{r}-\sqrt{r_{0}}-\sqrt{\log(1/\delta_{2})})^{2}}+\frac{2{r}\|{\mathcal{Z}}\|_{F}^{2}}{{\log(1/(1-\epsilon))}}.$
(8)

  3. 3.

(Low tubal-rank approximation with noisy sketches) Let $r_{1}$ be an integer
and define the error tensor $\mathcal{E}$ such that
$\mathcal{X}_{0}=\mathcal{X}_{1}+\mathcal{E}$ where $\mathcal{X}_{1}$ is the
best rank $r_{1}$ approximation of $\mathcal{X}_{0}$. If $r_{1}<r<n_{1}$ and
for all $k\in[n_{3}]$, $\hat{\tilde{\mathcal{Y}}}_{k}$ is of rank $r$, then
for any $\delta_{1},\delta_{2},\epsilon>0$ such that
$1\geq\delta_{2}>\exp{-(\sqrt{r}-\sqrt{r_{1}})^{2}}$ and $\epsilon<1$, with
probability at least $1-(\delta_{1}-\delta_{2}-\epsilon)n_{3}-2\delta_{2}$,

$\displaystyle\|\mathcal{X}-\mathcal{X}_{0}\|_{F}^{2}\leq
2\|\mathcal{E}\|^{2}_{F}\cdot$
$\displaystyle\left(\frac{4\sqrt{r(n_{1}-r)}(\sqrt{r}+\sqrt{n_{2}}+\sqrt{\log(1/\delta_{2})})}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})}+\frac{4\sqrt{r}\left(\sqrt{r}+\sqrt{n_{1}}+\sqrt{\log(1/\delta_{2})}\right)}{\sqrt{\log(1/(1-\epsilon))}}+1\right)^{2}$
$\displaystyle+\frac{{8r(n_{1}-r)}\|{\tilde{\mathcal{Z}}\|^{2}_{F}}}{{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})^{2}}+\frac{{8r}\|{\mathcal{Z}}\|_{F}^{2}}{{\log(1/(1-\epsilon))}}.$

The low tubal-rank approximation with noisy sketches error bound in Corollary
2 depends on $\|\mathcal{E}\|^{2}_{F}$ where
$\mathcal{E}=\mathcal{X}_{0}-\mathcal{X}_{1}$ and $\mathcal{X}_{1}$ is the
best tubal rank $r_{1}$ approximation of $\mathcal{X}_{0}$. It is well known
that the t-product satisfies a tensor Eckart-Young Theorem, i.e., the
$k$-truncated t-SVD produces the best tubal rank $k$ approximation of a given
matrix [23]. Thus the error term depends on the tail norms of the singular
tubes of each frontal slice of $\mathcal{X}_{0}$. In particular, we have that

$\displaystyle\|\mathcal{E}\|^{2}_{F}$
$\displaystyle=\sum_{i=r_{1}+1}^{\min(n_{1},n_{2})}\|S_{i,i,:}\|^{2},$ (9)

where $\mathcal{X}_{0}=\mathcal{U}\mathcal{S}\mathcal{V}^{T}$ is the t-SVD of
$\mathcal{X}_{0}$ as defined in Definition 1.

######  Remark 2.

In closely related work [31], the authors considered low-tubal-rank tensor
approximation from noiseless sketches, adapted and applied single-pass matrix
sketching algorithms from [38], while our setting here is to recover low-
tubal-rank tensors from noisy sketches. The results in [31] provided expected
error bounds in the noiseless case, while our work provides probabilistic
error bounds in the noisy case.

##  4 Proof of main results

Our proof for robust matrix recovery, presented in Theorem 2, derives an upper
bound for the difference between the output $X$ and the ground truth matrix
$X_{0}$. To accomplish this, the expression, $X-X_{0}$, is decomposed into two
summands. The first summand, which is responsible for the first error term in
(16), depends on $\tilde{Z}$ (and not on $Z$) and can be written as $PX_{0}$
for a projection matrix $P$. Next, we use the oblique projection matrix
expression in Lemma 3 to simplify the expression. We then control the error by
relating it to the smallest singular value of a truncated Haar unitary matrix.
Here we use the crucial fact that when $\tilde{Y}$ is full-rank,
$\textnormal{Im}(P)=\ker{S}$ is uniformly distributed on the Grassmannian of
all $(n-r)$-dimensional subspaces in $\mathbb{C}^{n}$. When $\tilde{Y}$ is not
full rank, $\textnormal{Im}(P)$ does not have that property, and our proof
technique cannot be directly applied. This part of the proof is summarized in
Lemma 1.

The second summand in the decomposition of $X-X_{0}$, which is responsible for
the second error term in (16), is due to the noise term $Z$, and is simpler to
handle. For this part, a lower bound on the smallest singular value of
Gaussian random matrices is utilized.

The distribution of the smallest singular value of truncated Haar unitary
matrices was explicitly calculated in [13, 5]. For a more general class of
random matrices (including Haar orthogonal matrices), such distribution was
derived in [13, 15, 4] in terms of generalized hypergeometric functions. By
using the corresponding tail probability bound [4, Corollary 3.4] for
truncated Haar orthogonal matrices, our analysis can be extended to real
Gaussian sketching matrices $S$ and $\tilde{S}$. Our proof techniques rely
heavily on the structure of Gaussian matrices with i.i.d. entries. For more
structured sketching matrices or sub-Gaussian sketching matrices, our current
approach can not be directly applied.

###  4.1 Proof of Theorem 1

###### Proof.

Let $X_{0}=U_{0}\Sigma_{0}V_{0}^{*}$ be the SVD of $X_{0}$ where $U_{0}$ is
$n_{1}\times r_{0}$, $V_{0}$ is $n_{2}\times r_{0}$ and $\Sigma_{0}$ is an
invertible, diagonal $r_{0}\times r_{0}$ matrix. Now we can write
$\overline{X}$ as

$\displaystyle\overline{X}$
$\displaystyle=U_{0}\Sigma_{0}V_{0}^{*}\tilde{S}^{*}(SU_{0}\Sigma_{0}V_{0}^{*}\tilde{S}^{*})^{\dagger}SU_{0}\Sigma_{0}V_{0}^{*}.$

Note that since $S,\tilde{S}$ are Gaussian matrices and since $U_{0}$ and
$V_{0}$ are orthonormal, $SU_{0}\in\mathbb{C}^{r\times r_{0}}$ and
$\tilde{S}V_{0}\in\mathbb{C}^{r\times r_{0}}$ have linearly independent
columns with probability 1 and by Property (1) in Lemma 2,

$(SU_{0})^{\dagger}SU_{0}=I,\quad
V_{0}^{*}\tilde{S}^{*}(V_{0}^{*}\tilde{S}^{*})^{\dagger}=I.$

So with probability 1,

$\displaystyle\overline{X}$
$\displaystyle=U_{0}\Sigma_{0}V_{0}^{*}\tilde{S}^{*}(V_{0}^{*}\tilde{S}^{*})^{\dagger}\Sigma_{0}^{-1}(SU_{0})^{\dagger}(SU_{0})\Sigma_{0}V_{0}^{*}=U_{0}\Sigma_{0}V_{0}^{*}=X_{0}.$

∎

###  4.2 Proof of Theorem 2

The double sketch algorithm outputs

$X=\tilde{Y}^{*}(S\tilde{Y}^{*})^{\dagger}(SX_{0}+Z)=\tilde{X}+\tilde{Y}^{*}(S\tilde{Y}^{*})^{\dagger}Z$

where

$\tilde{X}:=\tilde{Y}^{*}(S\tilde{Y}^{*})^{\dagger}SX_{0}.$ (10)

Let $A$ be such that

$A:=\tilde{Z}^{*}+X_{0}\tilde{S}^{*}\in\mathbb{C}^{n_{1}\times r}.$

Since $A=\tilde{Y}^{*}$ is of rank $r$ from the assumption in Theorem 2, we
denote the SVD of $A$ as

$\displaystyle A$ $\displaystyle=U_{A}\Sigma_{A}V_{A}^{*},$ $\displaystyle
U_{A}$ $\displaystyle\in\mathbb{C}^{n_{1}\times
r},\Sigma_{A}\in\mathbb{C}^{r\times r},V_{A}\in\mathbb{C}^{r\times r}.$ (11)

Furthermore, since $S$ is independent of $(\tilde{S},\tilde{Z})$, $SU_{A}$ is
invertible with probability $1$. Therefore,

$\displaystyle\tilde{Y}^{*}(S\tilde{Y}^{*})^{\dagger}$
$\displaystyle=A(SA)^{\dagger}$
$\displaystyle=(U_{A}\Sigma_{A}V_{A}^{*})(SU_{A}\Sigma_{A}V_{A}^{*})^{\dagger}$
$\displaystyle=U_{A}\Sigma_{A}V_{A}^{*}V_{A}\Sigma_{A}^{-1}(SU_{A})^{-1}$
$\displaystyle=U_{A}(SU_{A})^{-1},$ (12)

where in the third equality, we use the fact that $SU_{A}$ is invertible to
get

$\left(SU_{A}\Sigma_{A}V_{A}^{*}\right)^{\dagger}=V_{A}(SU_{A}\Sigma_{A})^{\dagger}=V_{A}\Sigma_{A}^{-1}(SU_{A})^{-1}.$

Using this notation, the output of (5) simplifies to

$X=\tilde{X}+\tilde{Y}^{*}(S\tilde{Y}^{*})^{\dagger}Z=\tilde{X}+U_{A}(SU_{A})^{-1}Z.$

Then

$X-X_{0}=\tilde{X}-X_{0}+\tilde{Y}^{*}(S\tilde{Y}^{*})^{-1}Z=(\tilde{X}-X_{0})+U_{A}(SU_{A})^{-1}Z,$

and since our goal is to bound the approximation error, we consider

$\displaystyle\vvvert X-X_{0}\vvvert$
$\displaystyle\leq\vvvert\tilde{X}-X_{0}\vvvert+\vvvert
U_{A}(SU_{A})^{-1}Z\vvvert.$

Lemma 1 allows us to bound the first term on the right-hand side of the
inequality above.

######  Lemma 1.

With probability at least $1-\delta_{1}-\delta_{2}$, where $\delta_{1}>0$, and
$1\geq\delta_{2}>\exp(-(\sqrt{r}-\sqrt{r_{0}})^{2})$, the matrix $\tilde{X}$
defined in equation (10) satisfies

$\displaystyle\vvvert\tilde{X}-X_{0}\vvvert\leq\frac{\sqrt{r(n_{1}-r)}\vvvert\tilde{Z}\vvvert}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{0}}-\sqrt{\log(1/\delta_{2})})}.$
(13)

###### Proof.

From (10) and (12),

$\displaystyle\tilde{X}-X_{0}$
$\displaystyle=\tilde{Y}^{*}(S\tilde{Y}^{*})^{\dagger}SX_{0}-X_{0}$
$\displaystyle=U_{A}(SU_{A})^{-1}SX_{0}-X_{0}$ (14)
$\displaystyle=-(I-U_{A}(SU_{A})^{-1}S)X_{0}$ $\displaystyle=-PX_{0},$

where we have set $P:=I-U_{A}(SU_{A})^{-1}S$. We observe that $P$ is a
projection, i.e., $P^{2}=P$, which satisfies

$\displaystyle\text{ker }P$ $\displaystyle=U_{A},$ $\displaystyle\text{Im }P$
$\displaystyle=\text{ker }S=:\tilde{V}\in\mathbb{C}^{n_{1}\times(n_{1}-r)}.$

Recall that $X_{0}=U_{0}\Sigma_{0}V_{0}^{*}$ is the SVD of $X_{0}$. From Lemma
3,

$\displaystyle PX_{0}$
$\displaystyle=\tilde{V}(U_{A,\perp}^{*}\tilde{V})^{-1}U_{A,\perp}^{*}X_{0}$
$\displaystyle=\tilde{V}(U_{A,\perp}^{*}\tilde{V})^{-1}U_{A,\perp}^{*}U_{0}\Sigma_{0}V_{0}^{*}.$

Then we can bound

$\displaystyle\vvvert\tilde{X}-X_{0}\vvvert$
$\displaystyle=\vvvert\tilde{V}(U_{A,\perp}^{*}\tilde{V})^{-1}U_{A,\perp}^{*}U_{0}\Sigma_{0}V_{0}^{*}\vvvert$
$\displaystyle\leq\|\tilde{V}\|_{2\rightarrow
2}\|(U_{A,\perp}^{*}\tilde{V})^{-1}\|_{2\rightarrow 2}\vvvert
U_{A,\perp}^{*}U_{0}\Sigma_{0}V_{0}^{*}\vvvert$ (15)
$\displaystyle=\|(U_{A,\perp}^{*}\tilde{V})^{-1}\|_{2\rightarrow 2}\vvvert
U_{A,\perp}^{*}U_{0}\Sigma_{0}V_{0}^{*}\vvvert.$ (16)

We focus our efforts on simplifying the second term of (16). Writing $A$ in
terms of the SVD of $X_{0}$, one obtains

$A=\tilde{Z}^{*}+X_{0}\tilde{S}^{*}=U_{A}\Sigma_{A}V_{A}^{*}=\tilde{Z}^{*}+U_{0}\Sigma_{0}V_{0}^{*}\tilde{S}^{*}.$

Rearranging terms yields

$U_{0}\Sigma_{0}V_{0}^{*}\tilde{S}^{*}=A-\tilde{Z}^{*}.$

Since $\tilde{S}V_{0}$ is of size $r\times r_{0}$, $\tilde{S}$ is complex
Gaussian and $V_{0}$ has orthonormal columns, then with probability $1$,
$\tilde{S}V_{0}$ has linearly independent columns. Thus
$V_{0}^{*}\tilde{S}^{*}$ has linearly independent rows and from Property (1)
in Lemma 2,

$(V_{0}^{*}\tilde{S}^{*})(V_{0}^{*}\tilde{S}^{*})^{\dagger}=I.$

This implies that

$U_{0}\Sigma_{0}=\left(A-\tilde{Z}^{*}\right)\left(V_{0}^{*}\tilde{S}^{*}\right)^{\dagger}.$

We note that

$\displaystyle\vvvert U_{A,\perp}^{*}U_{0}\Sigma_{0}V_{0}^{*}\vvvert$
$\displaystyle\leq\vvvert
U_{A,\perp}^{*}U_{0}\Sigma_{0}\vvvert\|V_{0}^{*}\|_{2\rightarrow 2}$
$\displaystyle=\vvvert U_{A,\perp}^{*}U_{0}\Sigma_{0}\vvvert$
$\displaystyle=\vvvert
U_{A,\perp}^{*}\left(A-\tilde{Z}^{*}\right)\left(V_{0}^{*}\tilde{S}^{*}\right)^{\dagger}\vvvert$
$\displaystyle=\vvvert
U_{A,\perp}^{*}\tilde{Z}^{*}\left(V_{0}^{*}\tilde{S}^{*}\right)^{\dagger}\vvvert$
$\displaystyle\leq\vvvert
U_{A,\perp}^{*}\tilde{Z}^{*}\vvvert\|\left(V_{0}^{*}\tilde{S}^{*}\right)^{\dagger}\|_{2\rightarrow
2}$ $\displaystyle=\frac{\vvvert
U_{A,\perp}^{*}\tilde{Z}^{*}\vvvert}{\sigma_{\min}\left(V_{0}^{*}\tilde{S}^{*}\right)}\leq\frac{\vvvert\tilde{Z}\vvvert}{\sigma_{\min}\left(V_{0}^{*}\tilde{S}^{*}\right)}.$

Using the fact that $\|(U_{A,\perp}^{*}\tilde{V})^{-1}\|_{2\rightarrow
2}=1/{\sigma_{\min}(U_{A,\perp}^{*}\tilde{V})}$ for the first term of (16), we
then obtain the following bound:

$\displaystyle\vvvert\tilde{X}-X_{0}\vvvert\leq\frac{\vvvert\tilde{Z}\vvvert}{\sigma_{\min}(U_{A,\perp}^{*}\tilde{V})\sigma_{\min}\left(\tilde{S}V_{0}\right)},$
(17)

which holds with probability 1.

We now derive a probabilistic bound from (17) using concentration inequalities
from random matrix theory. Since $\tilde{V}$ can be seen as the
$n_{1}\times(n_{1}-r)$ submatrix of a Haar unitary matrix
$(\tilde{V},\tilde{V}_{\perp})$ and $S$ is independent of
$(\tilde{Z},\tilde{S})$, by unitary invariance property,

$\displaystyle M:=\begin{bmatrix}U_{A,\perp}^{*}\\\
U_{A}^{*}\end{bmatrix}(\tilde{V},\tilde{V}_{\perp})$

is also a Haar unitary matrix, and $U_{A,\perp}^{*}\tilde{V}$ is exactly the
upper left $(n_{1}-r)\times(n_{1}-r)$ corner of $M$. We can apply Lemma 4 to
get

$\displaystyle\mathbb{P}\left(\|(U_{A,\perp}^{*}\tilde{V})^{-1}\|_{2\rightarrow
2}\leq\frac{\sqrt{r(n_{1}-r)}}{\sqrt{\delta}_{1}}\right)\geq 1-\delta_{1}$

for any $\delta_{1}>0$. Since $\tilde{S}V_{0}$ is distributed as an $r\times
r_{0}$ complex Gaussian random matrix, if $r>r_{0}$, by Lemma 5, for any
$\delta_{2}$ with
$1\geq\delta_{2}>\exp\left(-\left(\sqrt{r}-\sqrt{r_{0}}\right)^{2}\right)$,
with probability at least $1-\delta_{2}$,

$\displaystyle\sigma_{\min}(\tilde{S}V_{0})\geq\sqrt{r}-\sqrt{r_{0}}-\sqrt{\log(1/\delta_{2})}.$
(18)

Combining the two probability estimates, with probability at least
$1-\delta_{1}-\delta_{2}$,

$\displaystyle\vvvert\tilde{X}-X_{0}\vvvert\leq\frac{\sqrt{r(n_{1}-r)}\vvvert\tilde{Z}\vvvert}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{0}}-\sqrt{\log(1/\delta_{2})})}.$

∎

For the second term, by the independence between $S$ and $\tilde{Z}$, $SU_{A}$
is distributed as an $r\times r$ standard complex Gaussian matrix. Thus, it
follows from Lemma 6 that with probability at least $1-\epsilon$ for any
$0<\epsilon<1$,

$\sigma_{\min}(SU_{A})\geq\sqrt{\log(1/(1-\epsilon))}r^{-1/2}.$

Thus, with probability at least $1-\epsilon$,

$\displaystyle\vvvert U_{A}(SU_{A})^{-1}Z\vvvert$
$\displaystyle\leq\|U_{A}\|_{2\rightarrow 2}\|(SU_{A})^{-1}\|_{2\rightarrow
2}\vvvert Z\vvvert$ $\displaystyle=\frac{1}{\sigma_{\min}(SU_{A})}\vvvert
Z\vvvert$
$\displaystyle\leq\frac{\sqrt{r}}{\sqrt{\log(1/(1-\epsilon))}}\vvvert
Z\vvvert.$ (19)

Combining (13) and (19) yields the desired result. This finishes the proof of
Theorem 2.

###  4.3 Proof of Corollaries

###### Proof of Corollary 1.

Write $X_{0}=X_{1}+E$, where $X_{1}$ is the best rank-$r_{1}$ approximation to
$X_{0}$. We obtain

$\displaystyle Y$ $\displaystyle=SX_{0}+Z=SX_{1}+(SE+Z),$
$\displaystyle\tilde{Y}$
$\displaystyle=\tilde{S}X_{0}^{*}+\tilde{Z}=\tilde{S}X_{1}^{*}+(\tilde{S}E^{*}+\tilde{Z}).$

Applying Theorem 2 to $X_{1}$ and two error terms $(SE+Z)$ and
$\tilde{S}E^{*}+\tilde{Z}$, the following error bound holds with probability
at least $1-\delta_{1}-\delta_{2}-\epsilon$,

$\displaystyle\|X-X_{0}\|_{2\rightarrow 2}\leq\|X-X_{1}\|_{2\rightarrow
2}+\|E\|_{2\rightarrow 2}$
$\displaystyle\leq\frac{\sqrt{r(n_{1}-r)}\|\tilde{S}E^{*}+\tilde{Z}\|_{2\rightarrow
2}}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})}+\frac{\sqrt{r}\|SE+Z\|_{2\rightarrow
2}}{\sqrt{\log(1/(1-\epsilon))}}+\|E\|_{2\rightarrow 2}$
$\displaystyle\leq\left(\frac{\sqrt{r(n_{1}-r)}\|\tilde{S}\|_{2\rightarrow
2}}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})}+\frac{\sqrt{r}\|S\|_{2\rightarrow
2}}{\sqrt{\log(1/(1-\epsilon))}}+1\right)\|E\|_{2\rightarrow 2}$
$\displaystyle+\frac{\sqrt{r(n_{1}-r)}\|\tilde{Z}\|_{2\rightarrow
2}}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})}+\frac{\sqrt{r}\|Z\|_{2\rightarrow
2}}{\sqrt{\log(1/(1-\epsilon))}}$
$\displaystyle=\left(\frac{\sqrt{r(n_{1}-r)}\|\tilde{S}\|_{2\rightarrow
2}}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})}+\frac{\sqrt{r}\|S\|_{2\rightarrow
2}}{\sqrt{\log(1/(1-\epsilon))}}+1\right)\sigma_{r_{1}+1}\left(X_{0}\right)$
$\displaystyle+\frac{\sqrt{r(n_{1}-r)}\|\tilde{Z}\|_{2\rightarrow
2}}{\sqrt{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})}+\frac{\sqrt{r}\|Z\|_{2\rightarrow
2}}{\sqrt{\log(1/(1-\epsilon))}}.$

From concentration of the operator norm of $\tilde{S}$ and $S$ in Lemma 5,
with probability $1-2\delta_{2}$, we have that

$\|S\|_{2\rightarrow
2}\leq\sqrt{r}+\sqrt{n_{1}}+\sqrt{\log(1/\delta_{2})},\quad\|\tilde{S}\|_{2\rightarrow
2}\leq\sqrt{r}+\sqrt{n_{2}}+\sqrt{\log(1/\delta_{2})}.$

We then obtain the desired bound. ∎

###### Proof of Corollary 2.

Consider the double sketch tensor approach described in (7):

$\displaystyle\mathcal{Y}$
$\displaystyle=\mathcal{S}*\mathcal{X}_{0}+\mathcal{Z},\,\,\,\mathcal{\tilde{Y}}=\mathcal{\tilde{S}}*\mathcal{X}_{0}^{*}+\mathcal{\tilde{Z}}$

where $\mathcal{S}_{1}=S$, $\mathcal{\tilde{S}}_{1}=\tilde{S}$ and
$\mathcal{S}_{k}=\textbf{0}$, $\mathcal{\tilde{S}}_{k}=\textbf{0}$ for all
$k\in\\{2,\dots,n_{3}\\}$. With this construction, after performing a mode-3
Fourier transformation (2) on the tensors, the measurements
$\widehat{\mathcal{Y}}$ and $\widehat{\tilde{\mathcal{Y}}}$ can be decomposed
into $n_{3}$ low-rank matrix sketches:

$\displaystyle\widehat{\mathcal{Y}}_{i}$
$\displaystyle=\widehat{\mathcal{S}}_{i}\widehat{\mathcal{X}_{0}}_{i}+\widehat{\mathcal{Z}}_{i},\,\,\,\widehat{\mathcal{\tilde{Y}}}_{i}=\widehat{\mathcal{\tilde{S}}}_{i}\widehat{\mathcal{X}_{0}^{*}}_{i}+\widehat{\mathcal{\tilde{Z}}}_{i},\,\,\,\,i\in[n_{3}].$
(20)

Notice that it follows from the definition of mode-3 Fourier transformation in
(2), $\widehat{\mathcal{S}_{i}}$ and $\widehat{\tilde{\mathcal{S}_{i}}}$ are
complex Gaussian random matrices with independent entries by the unitary
invariance of the complex Gaussian distribution. Thus, the results of Theorem
1 and Theorem 2 can be applied to produce recovery guarantees for the double
sketch tensor approach. In particular, when $\mathcal{Z}=\mathbf{0}$,
$\tilde{\mathcal{Z}}=\mathbf{0}$, by Theorem 1, each of the $n_{3}$
transformed frontal slices $\mathcal{X}_{i}$ can be recovered with probability
1 from the sketches given in (20).

For general noise and when $r_{0}<r<n_{1}$, we start by noting that

$\displaystyle\|\mathcal{X}-\mathcal{X}_{0}\|_{F}^{2}$
$\displaystyle=\sum_{k=1}^{n_{3}}\|\mathcal{X}_{k}-(\mathcal{X}_{0})_{k}\|_{F}^{2}$
$\displaystyle\leq\sum_{k=1}^{n_{3}}\left(2\|\widehat{\mathcal{X}}_{k}-\widetilde{\widehat{\mathcal{X}}_{k}}\|_{F}^{2}+2\|\widetilde{\widehat{\mathcal{X}}_{k}}-(\widehat{\mathcal{X}}_{0})_{k}\|_{F}^{2}\right),$

where
$\widehat{\mathcal{X}}_{k},(\widehat{\mathcal{X}_{0}})_{k}\in\mathbb{C}^{n_{1}\times
n_{2}}$ are $k^{th}$ frontal slices of the mode-3 FFT of $\mathcal{X}$ and
$\mathcal{X}_{0}$, and $\widetilde{\widehat{\mathcal{X}}_{k}}$ is the
recovered frontal slice in the case where $\mathcal{Z}=0$. For the first term,
we can use the same argument as in the proof of Theorem 2 and invoke Lemma 1,
and for the second term, we apply (19). Thus, the approximation error is:

$\displaystyle\|\mathcal{X}-\mathcal{X}_{0}\|_{F}^{2}$
$\displaystyle\leq\sum_{k=1}^{n_{3}}\left(2\|\widehat{\mathcal{X}}_{k}-\widetilde{\widehat{\mathcal{X}}_{k}}\|_{F}^{2}+2\|\widetilde{\widehat{\mathcal{X}}_{k}}-(\widehat{\mathcal{X}}_{0})_{k}\|_{F}^{2}\right)$
$\displaystyle\leq\sum_{k=1}^{n_{3}}\left(\frac{2{r(n_{1}-r)}\|\widehat{\tilde{\mathcal{Z}_{k}}}\|_{F}^{2}}{{\delta_{1}}(\sqrt{r}-\sqrt{r_{0}}-\sqrt{\log(1/\delta_{2})})^{2}}+\frac{{2r}\|\widehat{\mathcal{Z}_{k}}\|_{F}^{2}}{{\log(1/(1-\epsilon))}}\right)$
$\displaystyle=\frac{2{r(n_{1}-r)}\|{\tilde{\mathcal{Z}}}\|_{F}^{2}}{{\delta_{1}}(\sqrt{r}-\sqrt{r_{0}}-\sqrt{\log(1/\delta_{2})})^{2}}+\frac{2{r}\|{\mathcal{Z}}\|_{F}^{2}}{{\log(1/(1-\epsilon))}}.$
(21)

We now prove the low-tubal-rank approximation case, which follows the proof of
Corollary 1. Let $r_{1}$ be an integer and define the error tensor
$\mathcal{E}$ such that $\mathcal{X}_{0}=\mathcal{X}_{1}+\mathcal{E}$ where
$\mathcal{X}_{1}$ is the best rank $r_{1}$ approximation of $\mathcal{X}_{0}$.
Thus, the sketches can be written as

$\displaystyle\mathcal{Y}$
$\displaystyle=\mathcal{S}*\mathcal{X}_{0}+\mathcal{Z}=\mathcal{S}*\mathcal{X}_{1}+(\mathcal{S}*\mathcal{E}+\mathcal{Z})$
$\displaystyle\mathcal{\widetilde{Y}}$
$\displaystyle=\mathcal{\widetilde{S}}*\mathcal{X}_{0}^{*}+\mathcal{\widetilde{Z}}=\mathcal{\widetilde{S}}*\mathcal{X}_{1}^{*}+(\mathcal{\widetilde{S}}*\mathcal{E}^{*}+\mathcal{\widetilde{Z}}),$

and thus we can interpret $\mathcal{Y}$ and $\mathcal{\tilde{Y}}$ as sketches
of $\mathcal{X}_{1}$. If $r_{1}<r<n_{1}$ and for all $k\in[n_{3}]$,
$\hat{\tilde{\mathcal{Y}}}_{k}$ is of rank $r$, we have that

$\displaystyle\|\mathcal{X}-\mathcal{X}_{0}\|^{2}_{F}$ $\displaystyle\leq
2\|\mathcal{X}-\mathcal{X}_{1}\|^{2}_{F}+2\|\mathcal{E}\|^{2}_{F}$
$\displaystyle\leq\frac{4{r(n_{1}-r)}\|\mathcal{\widetilde{S}}*\mathcal{E}^{*}+\mathcal{\widetilde{Z}}\|_{F}^{2}}{{\delta_{1}}(\sqrt{r}-\sqrt{r_{1}}-\sqrt{\log(1/\delta_{2})})^{2}}+\frac{4{r}\|\mathcal{S}*\mathcal{E}+\mathcal{Z}\|_{F}^{2}}{{\log(1/(1-\epsilon))}}+2\|\mathcal{E}\|^{2}_{F},$
(22)

with probability at least $1-(\delta_{1}-\delta_{2}-\epsilon)n_{3}$, where the
second inequality uses (21). We can further decompose the term:

$\displaystyle\|\mathcal{S}*\mathcal{E}+\mathcal{Z}\|_{F}^{2}$
$\displaystyle\leq
2\|\mathcal{S}*\mathcal{E}\|_{F}^{2}+2\|\mathcal{Z}\|_{F}^{2}$
$\displaystyle=2\sum_{k=1}^{n_{3}}\|\widehat{\mathcal{S}}_{k}\widehat{\mathcal{E}}_{k}\|^{2}_{F}+2\|\mathcal{Z}\|_{F}^{2}$
$\displaystyle\leq
2\sum_{k=1}^{n_{3}}\|\widehat{\mathcal{S}}_{k}\|^{2}_{2\rightarrow
2}\|\widehat{\mathcal{E}}_{k}\|^{2}_{F}+2\|\mathcal{Z}\|_{F}^{2}$
$\displaystyle=2\|S\|^{2}_{2\rightarrow
2}\|{\mathcal{E}}\|^{2}_{F}+2\|\mathcal{Z}\|_{F}^{2}~{},$

where we use the fact that $\widehat{\mathcal{S}}_{k}=S$ for all $k$ in the
last equality. By Lemma 5, we have that with probability at least
$1-\delta_{2}$,

$\displaystyle\|\mathcal{S}*\mathcal{E}+\mathcal{Z}\|_{F}^{2}$
$\displaystyle\leq
2\left(\sqrt{r}+\sqrt{n_{1}}+\sqrt{\log(1/\delta_{2})}\right)\|{\mathcal{E}}\|^{2}_{F}+2\|\mathcal{Z}\|_{F}^{2}.$

Similarly, for
$\|\mathcal{\widetilde{S}}\mathcal{E}^{*}+\mathcal{\widetilde{Z}}\|_{F}^{2}$,
we have that:

$\displaystyle\|\mathcal{\widetilde{S}}\mathcal{E}^{*}+\mathcal{\widetilde{Z}}\|_{F}^{2}$
$\displaystyle\leq
2\left(\sqrt{r}+\sqrt{n_{2}}+\sqrt{\log(1/\delta_{2})}\right)\|\mathcal{E}\|^{2}_{F}+2\|\widetilde{\mathcal{Z}}\|^{2}_{F},$

with probability $1-\delta_{2}$. Finally, using the upper bounds on
$\|\mathcal{\widetilde{S}}\mathcal{E}^{*}+\mathcal{\widetilde{Z}}\|_{F}^{2}$
and $\|\mathcal{S}\mathcal{E}+\mathcal{Z}\|_{F}^{2}$ in (22) obtains the
desired result. ∎

##  5 Experiments

We complement our theoretical findings with numerical experiments examining
the empirical behavior of the approximation error using the double sketch
algorithm for recovering low-rank matrices from noise sketches111Code is
available at: https://github.com/anna-math/doublesketch. In these experiments,
we fix the dimension of the rank $r_{0}$ matrix
$X_{0}\in\mathbb{R}^{n_{1}\times n_{2}}$ to be $n=n_{1}=n_{2}=100$,
initialized by multiplying two i.i.d. complex Gaussian matrices of size
$100\times r_{0}$ and $r_{0}\times 100$. (The singular value distribution of
the product of two rectangular random Gaussian matrices can be found in [8].
Notably, this is not the Marchenko-Pastur law [3]. Empirical approximation
error represents the median error over 50 trials where the sketching matrices
change with each trial. The entries of the error matrices (or tensors) $Z$ and
$\tilde{Z}$ ($\mathcal{Z}$ and $\tilde{\mathcal{Z}}$) are drawn i.i.d. from a
standard Gaussian distribution and then normalized to achieve the desired
Frobenius norm. For these numerical experiments, all norms and errors are
measured in the Frobenius norm.

Figure 1 demonstrates the performance of the double sketch algorithm for
recovering a low-rank matrix when $r_{0}=10$ and
$r\in\\{r_{0}+1,2r_{0},n-1\\}$ (left to right). For each tuple $(r,r_{0})$, we
plot the approximation error while varying the noise levels
$\epsilon_{1}=\|Z\|_{F}$ and $\epsilon_{2}=\|\tilde{Z}\|_{F}$. Note that in
Figure 1 as $r$ increases and approaches $n_{1}$, the impact of $\epsilon_{2}$
in the approximation error becomes weaker. This can be seen best in the far
right plot where the choice of $\epsilon_{2}$ seems to have no influence on
the approximation error of the double sketch algorithm. This is expected as
the error bound represented in Theorem 2 depends on $n-r$ and if we pick
$r=n$, the error bound becomes independent of $\epsilon_{2}$.

Figure 1: Numerical results demonstrating the empirical performance measured
by the approximation error $\|X-X_{0}\|_{F}$ of the double sketching algorithm
for recovering a $100\times 100$ rank $r_{0}=10$ matrix while varying
$\|{{Z}}\|_{F}=\epsilon_{1}$ and $\|{\tilde{Z}}\|_{F}=\epsilon_{2}$ when
(left) $r=r_{0}+1$ (center) $r=2r_{0}$ (right) $r=n-1$. Approximation error
reported is the median after 50 trials.

In our next experiment, we consider the recovery of low-tubal-rank tensors.
Figure 2 presents the empirical approximation errors obtained when using the
double sketch Algorithm 7 to recover low-tubal-rank tensors. To initialize the
underlying low-rank tensor, we multiply two tensors of dimension $n\times
r_{0}\times n_{3}$ and $r_{0}\times n\times n_{3}$ such that
$r_{0}<\min\\{n,n_{3}\\}$ and we set $n=100$ in this experiment. Thus, the
tubal rank of the tensor is $r_{0}$ and $\mathcal{X}_{0}\in\mathbb{R}^{n\times
n\times n_{3}}$. Similar to the experiment presented for the matrix case, the
approximation error is the median error from 50 trials where new sketching
matrices are drawn for each trial. In the tensor case, as we vary
$\epsilon_{1}$ and $\epsilon_{2}$, the behavior of the approximation error is
similar to the matrix case. In particular, as $r$ approaches $n$, the error
bound becomes independent of $\epsilon_{2}$. In other words, the trend for the
approximation error does not change as $\epsilon_{2}$ varies, as shown in the
far right plot of Figure 2.

Figure 2: The Frobenius norm error $\|\mathcal{X}-\mathcal{X}_{0}\|$ of the
double sketching algorithm for recovering a $100\times 100\times n_{3}$ tensor
of tubal-rank $r_{0}=10$ varying $\|\mathcal{\tilde{Z}}\|_{F}=\epsilon_{1}$
and $\|\mathcal{{Z}}\|_{F}=\epsilon_{2}$ for: (left column) $r=r_{0}+1$
(center column) $r=2r_{0}$ (right column) $r=n-1$.

Figure 3 presents the performance of Algorithm 2 as $n_{3}$ increases for
difference choices of $r$. For this experiment, $\mathcal{X}_{0}$ is a tubal-
rank $r_{0}=10$ tensor of dimension $100\times 100\times n_{3}$ with unit
Frobenius norm. The noise level is set to $\|\mathcal{\tilde{Z}}\|_{F}=0.01$
and $\|\mathcal{{Z}}\|_{F}=0.01$. Here we make two observations. First, we
note that as $r$ approaches $n=100$, the approximation error gets smaller.
This is expected as the first term of the error bound (8) depending on
$\|\mathcal{\tilde{Z}}\|_{F}$ decreases as $r$ approaches $n$. Second, as
$n_{3}$ increases, we see that the approximation error slowly increases with
$n_{3}$, which is expected since our probability bound gets worse as $n_{3}$
increases.

Figure 3: Frobenius norm error of the double sketching algorithm for
recovering a $100\times 100\times n_{3}$ while varying $n_{3}$. The noise
level is set to $\|\mathcal{\tilde{Z}}\|_{F}=0.01$ and
$\|\mathcal{{Z}}\|_{F}=0.01$.

In our final experiment, we demonstrate the performance of Algorithm 2 on CT
scan data. This data set contains CT scan slices of the C1 vertebrae. The
image datasets used in this experiment were from the Laboratory of Human
Anatomy and Embryology, University of Brussels (ULB), Belgium222Bone and joint
CT-Scan Data https://isbweb.org/data/vsj/. Here the ground-truth tensor is
given by $\mathcal{X}_{0}\in\mathbb{R}^{512\times 512\times 47}$. We pick
$r=200$ and compare the application of Algorithm 2 to the application of
Algorithm 1 naively to each slice of the tensor (which we refer to as slice-
wise matrix recovery). Figure 4 presents a visualization of the recovery of
the 20th slice of the tensor. The first image is the ground truth slice, the
second image is the 20th frame treated as a matrix and using Algorithm 1 for
recovery, and the third image is the 20th frame recovered from the double
sketched tensor using Algorithm 2. The approximation error
$\|\mathcal{X}-\mathcal{X}_{0}\|_{F}$ using Algorithm 2 is $0.3825$ whereas
the approximation error applying Algorithm 1 is $0.8541$. Algorithm 2 has the
added benefit of only requiring a single sketching matrix
$S\in\mathbb{R}^{512\times 200}$ whereas naively applying Algorithm 1 would
require $n_{3}$ sketching matrices to be stored. Of course, it would be
natural to ask whether one can simply use the same sketching matrix in
Algorithm 1. Doing so attained a worse approximation error of $1.6181$.

Figure 4: Performance of Algorithm 1 applied frontal slice-wise to a tensor
$\mathcal{X}_{0}\in\mathbb{R}^{512\times 512\times 47}$ containing CT scan
data of the C1 vertebrae. Here, we choose $r=200$ and
$\epsilon_{1}=\epsilon_{2}=0.01$. This figure shows the 20th frontal slice for
(left) ground truth tensor (center) recovered frontal slice after applying
Algorithm 1 to each frontal slice (right) recovered frontal slice after
applying Algorithm 2 to recover the ground truth tensor.

##  6 Conclusion

In this paper, we study the problem of recovering low-rank matrices and low-
tubal-rank tensors from two sketched matrices (or tensors), which are
perturbed by additive noise matrices, which we denote by $Z$, respectively
$\tilde{Z}$. This single-pass method has been proposed, but not theoretically
analyzed, in [16]. Our main results show that, without making assumptions on
the distribution of the noise matrices $Z$ and $\tilde{Z}$, one can
reconstruct the low-rank matrix using the method discussed in [16] and obtain
theoretical guarantees when the sketching matrices are assumed to be Gaussian
matrices. Furthermore, we show that this analysis is also applicable to low-
tubal-rank tensor recovery. Numerical experiments corroborate our theoretical
findings.

While in this paper we have used the original sketching method considered in
[16], recent years have seen the development of sketching methods for low-rank
approximation methods, which are, in general, more accurate. For example, the
method in [39] used an additional third matrix sketch, which is called core
sketch, to improve the accuracy. It would also be interesting to analyze the
noise robustness of these more recent methods and to compare the results with
the method analyzed in this paper.

Another possible extension is to analyze the truncated version of the noisy
sketching method in [16, 46] for low-rank approximation, which utilizes a
truncated SVD on the sketches $Y$ and $\tilde{Y}$ before performing the
recovery of $X_{0}$ (5).

Beyond sketching for matrices and low tubal-rank tensors, it would be
interesting to see whether we can analyze the reconstruction of low Tucker-
rank tensors from noisy sketches. A low Tucker-rank approximation of an
order-$N$ Tucker-rank tensor using linear sketches was considered in prior
work [36], where the authors used $N$ random sketching tensors to obtain $N$
factor sketches, and another $N$ random sketching matrices to obtain a core
sketch. Based on those sketches, a single-pass low Tucker-rank approximation
method was studied. In the matrix case ($N=2$), their algorithm corresponds to
the one in [39]. Since the Tucker decomposition depends on tensor unfolding
operations, and the sketching algorithms are applied to the unfolded tensors,
it would be interesting to extend our noisy double sketching analysis to their
setting (see Algorithm 4.2 in [36]).

It will also be interesting to see whether our analysis can be extended to the
scenario where more structured sketching matrices such as subsampled
randomized Hadamard transforms [7] or subsampled randomized Fourier transforms
(SRFT) are used. These more structured sketching matrices have the advantage
of reducing the computational burden, and the memory footprint compared to
Gaussian sketching matrices (see, e.g., [38] for sketching algorithms for low-
rank matrix approximation using SRFT matrices). Note that while our
deterministic formula for the difference between the output of the sketching
algorithm and the ground truth matrix also holds in this setting, it is not
known to us whether some of the involved concentration inequalities still hold
in this scenario. We believe that this is an interesting avenue for future
research.

##  Appendix A Technical lemmas

In this section, we collect a few technical lemmas, which were needed in our
proofs.

######  Lemma 2 (Properties of the pseudo-inverse, see, e.g., [17]).

The following properties hold for the pseudo-inverse.

  1. 1.

If $A$ has linearly independent columns, then $A^{\dagger}A=I$. If $B$ has
linearly independent rows, then $BB^{\dagger}=I$.

  2. 2.

$(AB)^{\dagger}=B^{\dagger}A^{\dagger}$ if $A$ has orthonormal columns or $B$
has orthonormal rows.

  3. 3.

$(AB)^{\dagger}=B^{\dagger}A^{\dagger}$ if $A$ has linearly independent
columns.

  4. 4.

If $A$ has orthonormal columns or orthonormal rows, then $A^{\dagger}=A^{*}$.

######  Lemma 3 (Oblique projection matrix, Theorem 2.1 in [20], see also
Equation (7.10.40) in [27]).

Let $P\in\mathbb{C}^{n\times n}$ be a projection matrix, i.e., $P^{2}=P$. We
note that $P$ is uniquely characterized by the subspaces $V_{1}$ and $V_{2}$
given by

$\displaystyle V_{1}$ $\displaystyle:=\textnormal{Im }P,$ $\displaystyle
V_{2}$ $\displaystyle=\textnormal{Ker }P.$

and $\text{dim }V_{1}+\text{dim }V_{2}=n$. By some abuse of notation, we
denote by $V_{1}$ a matrix with orthonormal columns, whose column span is
equal to the subspace $V_{1}$. Analogously, we denote by $V_{2}$ a matrix with
orthonormal columns, whose column span is equal to the subspace $V_{2}$. It
holds that

$\displaystyle P=V_{1}(V_{2,\perp}^{*}V_{1})^{-1}V_{2,\perp}^{*}.$

######  Lemma 4 (Smallest singular value of a truncated Haar unitary matrix,
Proposition C.3 in [5]).

Let $A$ be the upper left $(n-r)\times(n-r)$ corner of a Haar unitary matrix
$U$. Then for any $\delta>0$

$\displaystyle\mathbb{P}\left(\sigma_{\min}(A)\geq\frac{\sqrt{\delta}}{\sqrt{r(n-r)}}\right)\geq
1-\delta.$

######  Lemma 5 (Extreme singular value of Gaussian random matrices, Corollary
5.35 in [41]).

Let $A$ be an $m\times n$ random matrix with independent standard complex
Gaussian entries and $m>n$, then for any $\delta$ with $0<\delta\leq 1$, with
probability at least $1-\delta$,

$\sigma_{\min}(A)\geq\sqrt{m}-\sqrt{n}-\sqrt{\log(1/\delta)}.$

Similarly, with probability at least $1-\delta$,

$\sigma_{\max}(A)\leq\sqrt{m}+\sqrt{n}+\sqrt{\log(1/\delta)}.$

###### Proof.

This lemma was proved in [41] for real Gaussian matrices. We include a proof
for complex Gaussian matrices for completeness. We have

$\displaystyle\sigma_{\min}(A)=\min_{u\in S^{n-1}}\max_{v\in S^{m-1}}\langle
Au,v\rangle,$

where $S^{n-1}$ is the unit sphere in $\mathbb{C}^{n}$. Let $X_{u,v}=\langle
Au,v\rangle$ and $Y_{u,v}=\langle g,u\rangle+\langle h,v\rangle$, where
$g\in\mathbb{C}^{n}$, $h\in\mathbb{C}^{m}$ are independent standard complex
Gaussian random vectors. We find for any $x,w\in S^{n-1}$ and any $y,z\in
S^{m-1}$,

$\displaystyle\mathbb{E}|X_{u,v}-X_{w,z}|^{2}=\|uv^{*}-wz^{*}\|_{F}^{2}\leq\|u-w\|_{2}^{2}+\|v-z\|_{2}^{2}=\mathbb{E}|Y_{u,v}-Y_{w,z}|^{2}.$

From Gordon’s inequality ([42, Exercise 7.2.14]) and [42, Exercise 7.3.4],

$\mathbb{E}\sigma_{\min}(A)\geq\mathbb{E}\|h\|_{2}-\mathbb{E}\|g\|_{2}.$

We also have $\|h\|_{2}^{2}=\sum_{i=1}^{m}|h_{i}|^{2}=\frac{1}{2}W$, where
$W\sim\chi^{2}(2m)$. So

$\displaystyle\mathbb{E}\|h\|_{2}-\sqrt{m}=\frac{1}{\sqrt{2}}\left(\mathbb{E}\sqrt{W}-\sqrt{2m}\right).$

Note that $f(n)=\mathbb{E}\|g\|_{2}-\sqrt{n}$ is an increasing function in
$n$, where $g\sim N(0,I_{n})$, and $\|g\|_{2}^{2}\sim\chi^{2}(n)$. So for
$m\geq n$,

$\mathbb{E}\|h\|_{2}-\sqrt{m}\geq\mathbb{E}\|g\|_{2}-\sqrt{n},$

which gives $\mathbb{E}\sigma_{\min}(A)\geq\sqrt{m}-\sqrt{n}$. Since
$\sigma_{\min}(A)$ is a Lipschitz function of $2mn$ independent real Gaussian
random variables with mean zero and variance $1/2$, the probability estimate
follows from Gaussian concentration of Lipschitz functions (see [41,
Proposition 5.34]). The largest singular value $\sigma_{\max}(A)$ can be
estimated in a similar way using the Sudakov-Fernique inequality see, e.g.,
[42, Theorem 7.2.11]. ∎

######  Lemma 6 (Smallest eigenvalue of a square Gaussian random matrix, [14]
and Theorem 1.1 in [37]).

Let $A$ be an $n\times n$ complex standard Gaussian random matrix. Then, for
any $\epsilon>0$,

$\displaystyle\mathbb{P}\left(\sigma_{\min}(A)\geq\epsilon
n^{-1/2}\right)=e^{-\epsilon^{2}}.$

## References

  * [1] Raman Arora and Jalaj Upadhyay.  Differentially private robust low-rank approximation.  Advances in neural information processing systems, 31, 2018.
  * [2] Martin Azizyan, Akshay Krishnamurthy, and Aarti Singh.  Extreme compressive sampling for covariance estimation.  IEEE Trans. Inf. Theory, 64(12):7613–7635, 2018.
  * [3] Zhidong Bai and Jack W Silverstein.  Spectral analysis of large dimensional random matrices, volume 20.  Springer, 2010.
  * [4] Grey Ballard, James Demmel, Ioana Dumitriu, and Alexander Rusciano.  A generalized randomized rank-revealing factorization.  arXiv preprint arXiv:1909.06524, 2019.
  * [5] Jess Banks, Jorge Garza-Vargas, Archit Kulkarni, and Nikhil Srivastava.  Pseudospectral shattering, the sign function, and diagonalization in nearly matrix multiplication time.  Found. Comput. Math., pages 1–89, 2022.
  * [6] Ayan Kumar Bhunia, Subhadeep Koley, Abdullah Faiz Ur Rahman Khilji, Aneeshan Sain, Pinaki Nath Chowdhury, Tao Xiang, and Yi-Zhe Song.  Sketching without worrying: Noise-tolerant sketch-based image retrieval.  In Proceedings of the IEEE/CVF Conference on Computer Vision and Pattern Recognition, pages 999–1008, 2022.
  * [7] Christos Boutsidis and Alex Gittens.  Improved matrix algorithms via the subsampled randomized Hadamard transform.  SIAM J. Matrix Anal. Appl., 34(3):1301–1340, 2013.
  * [8] Zdzisław Burda, Andrzej Jarosz, Giacomo Livan, Maciej A Nowak, and Artur Swiech.  Eigenvalues and singular values of products of rectangular gaussian random matrices.  Phys. Rev. E, 82(6):061114, 2010.
  * [9] Venkat Chandrasekaran, Benjamin Recht, Pablo A Parrilo, and Alan S Willsky.  The convex geometry of linear inverse problems.  Found. Comput. Math., 12(6):805–849, 2012.
  * [10] Yuxin Chen, Yuejie Chi, Jianqing Fan, and Cong Ma.  Spectral methods for data science: A statistical perspective.  Foundations and Trends® in Machine Learning, 14(5):566–806, 2021.
  * [11] Kenneth L Clarkson and David P Woodruff.  Numerical linear algebra in the streaming model.  In Proceedings of the forty-first annual ACM symposium on Theory of computing, pages 205–214, 2009.
  * [12] Chandler Davis and William Morton Kahan.  The rotation of eigenvectors by a perturbation. iii.  SIAM J. Numer. Anal., 7(1):1–46, 1970.
  * [13] Ioana Dumitriu.  Smallest eigenvalue distributions for two classes of $\beta$-jacobi ensembles.  J. Math. Phys., 53(10):103301, 2012.
  * [14] Alan Edelman.  Eigenvalues and condition numbers of random matrices.  SIAM J. Matrix Anal. Appl., 9(4):543–560, 1988.
  * [15] Alan Edelman and Brian D Sutton.  The beta-Jacobi matrix model, the CS decomposition, and generalized singular value problems.  Found. Comput. Math., 8(2):259–285, 2008.
  * [16] Maryam Fazel, E Candes, Benjamin Recht, and P Parrilo.  Compressed sensing and robust recovery of low rank matrices.  In 2008 42nd Asilomar Conference on Signals, Systems and Computers, pages 1043–1047. IEEE, 2008.
  * [17] G.H. Golub and C.F. Van Loan.  Matrix Computations.  Johns Hopkins Studies in the Mathematical Sciences. Johns Hopkins University Press, 2013.
  * [18] Nathan Halko, Per-Gunnar Martinsson, Yoel Shkolnisky, and Mark Tygert.  An algorithm for the principal component analysis of large data sets.  SIAM J. Sci. Comput., 33(5):2580–2594, 2011.
  * [19] Nathan Halko, Per-Gunnar Martinsson, and Joel A Tropp.  Finding structure with randomness: Probabilistic algorithms for constructing approximate matrix decompositions.  SIAM Rev., 53(2):217–288, 2011.
  * [20] Per Christian Hansen.  Oblique projections, pseudoinverses, and standard-form transformations.  Technical University of Denmark, Tech. Rep, 2004.
  * [21] Botao Hao, Anru R Zhang, and Guang Cheng.  Sparse and low-rank tensor estimation via cubic sketchings.  In International Conference on Artificial Intelligence and Statistics, pages 1319–1330. PMLR, 2020.
  * [22] Eric Kernfeld, Misha Kilmer, and Shuchin Aeron.  Tensor–tensor products with invertible linear transforms.  Linear Algebra Appl., 485:545–570, 2015.
  * [23] Misha E Kilmer, Karen Braman, Ning Hao, and Randy C Hoover.  Third-order tensors as operators on matrices: A theoretical and computational framework with applications in imaging.  SIAM J. Matrix Anal. Appl., 34(1):148–172, 2013.
  * [24] Misha E Kilmer and Carla D Martin.  Factorization strategies for third-order tensors.  Linear Algebra Appl., 435(3):641–658, 2011.
  * [25] Canyi Lu, Jiashi Feng, Yudong Chen, Wei Liu, Zhouchen Lin, and Shuicheng Yan.  Tensor robust principal component analysis with a new tensor nuclear norm.  IEEE transactions on pattern analysis and machine intelligence, 42(4):925–938, 2019.
  * [26] Canyi Lu, Jiashi Feng, Zhouchen Lin, and Shuicheng Yan.  Exact low tubal rank tensor recovery from gaussian measurements.  In Proceedings of the 27th International Joint Conference on Artificial Intelligence, pages 2504–2510, 2018.
  * [27] Carl D. Meyer.  Matrix analysis and applied linear algebra.  Philadelphia, PA: Society for Industrial and Applied Mathematics (SIAM), 2000.
  * [28] Shanmugavelayutham Muthukrishnan.  Data streams: Algorithms and applications.  Foundations and Trends® in Theoretical Computer Science, 1(2):117–236, 2005.
  * [29] Sahand Negahban and Martin J. Wainwright.  Estimation of (near) low-rank matrices with noise and high-dimensional scaling.  Ann. Stat., 39(2):1069–1097, 2011.
  * [30] Sean O’Rourke, Van Vu, and Ke Wang.  Random perturbation of low rank matrices: Improving classical bounds.  Linear Algebra Appl., 540:26–59, 2018.
  * [31] Liqun Qi and Gaohang Yu.  T-singular values and t-sketching for third order tensors.  arXiv preprint arXiv:2103.00976, 2021.
  * [32] Benjamin Recht, Maryam Fazel, and Pablo A. Parrilo.  Guaranteed minimum-rank solutions of linear matrix equations via nuclear norm minimization.  SIAM Rev., 52(3):471–501, 2010.
  * [33] Mark Rudelson and Roman Vershynin.  Non-asymptotic theory of random matrices: extreme singular values.  In Proceedings of the International Congress of Mathematicians 2010 (ICM 2010) (In 4 Volumes) Vol. I: Plenary Lectures and Ceremonies Vols. II–IV: Invited Lectures, pages 1576–1602. World Scientific, 2010.
  * [34] Rakshith Sharma Srinivasa, Kiryung Lee, Marius Junge, and Justin Romberg.  Decentralized sketching of low rank matrices.  In H. Wallach, H. Larochelle, A. Beygelzimer, F. d'Alché-Buc, E. Fox, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 32. Curran Associates, Inc., 2019.
  * [35] Rakshith Sharma Srinivasa, Kiryung Lee, Marius Junge, and Justin Romberg.  Decentralized sketching of low rank matrices.  Advances in Neural Information Processing Systems, 32, 2019.
  * [36] Yiming Sun, Yang Guo, Charlene Luo, Joel Tropp, and Madeleine Udell.  Low-rank tucker approximation of a tensor from streaming data.  SIAM Journal on Mathematics of Data Science, 2(4):1123–1150, 2020\.
  * [37] Terence Tao and Van Vu.  Random matrices: The distribution of the smallest singular values.  Geom. Funct. Anal., 20(1):260–297, 2010.
  * [38] Joel A. Tropp, Alp Yurtsever, Madeleine Udell, and Volkan Cevher.  Practical sketching algorithms for low-rank matrix approximation.  SIAM J. Matrix Anal. Appl., 38(4):1454–1485, 2017.
  * [39] Joel A. Tropp, Alp Yurtsever, Madeleine Udell, and Volkan Cevher.  Streaming low-rank matrix approximation with an application to scientific simulation.  SIAM J. Sci. Comput., 41(4):a2430–a2463, 2019.
  * [40] Jalaj Upadhyay.  The price of privacy for low-rank factorization.  Advances in Neural Information Processing Systems, 31, 2018.
  * [41] Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices, page 210–268.  Cambridge University Press, 2012.
  * [42] Roman Vershynin.  High-dimensional probability: An introduction with applications in data science, volume 47.  Cambridge university press, 2018.
  * [43] Hailin Wang, Feng Zhang, Jianjun Wang, Tingwen Huang, Jianwen Huang, and Xinling Liu.  Generalized nonconvex approach for low-tubal-rank tensor recovery.  IEEE Trans. Neural Netw. Learn. Syst., 2021.
  * [44] Per-Åke Wedin.  Perturbation bounds in connection with singular value decomposition.  BIT, 12(1):99–111, 1972.
  * [45] Bernard Widrow and István Kollár.  Quantization noise: roundoff error in digital computation, signal processing, control, and communications.  Cambridge University Press, 2008.
  * [46] Franco Woolfe, Edo Liberty, Vladimir Rokhlin, and Mark Tygert.  A fast randomized algorithm for the approximation of matrices.  Appl. Comput. Harmon. Anal., 25(3):335–366, 2008.
  * [47] Qianxin Yi, Chenhao Wang, Kaidong Wang, and Yao Wang.  Effective streaming low-tubal-rank tensor approximation via frequent directions.  IEEE Transactions on Neural Networks and Learning Systems, 2022\.
  * [48] Feng Zhang, Wendong Wang, Jingyao Hou, Jianjun Wang, and Jianwen Huang.  Tensor restricted isometry property analysis for a large class of random measurement ensembles.  Science China Information Sciences, 64(1):1–3, 2021.
  * [49] Feng Zhang, Wendong Wang, Jianwen Huang, Jianjun Wang, and Yao Wang.  Rip-based performance guarantee for low-tubal-rank tensor recovery.  J. Comput. Appl. Math., 374:112767, 2020.
  * [50] Zemin Zhang and Shuchin Aeron.  Exact tensor completion using t-SVD.  IEEE Trans. Signal Process., 65(6):1511–1526, 2016.
  * [51] Tianyi Zhou and Dacheng Tao.  Bilateral random projections.  In 2012 IEEE International Symposium on Information Theory Proceedings, pages 1286–1290. IEEE, 2012.