# The Product of Gaussian Matrices is Close to Gaussian

Yi Li111Supported in part by Singapore Ministry of Education (AcRF) Tier 2
grant MOE2018-T2-1-013
Division of Mathematical Sciences
Nanyang Technological University
<EMAIL_ADDRESS>David P. Woodruff222Supported in part by Office of Naval
Research (ONR) grant N00014-18-1-256 and a Simons Investigator Award
Department of Computer Science
Carnegie Mellon University
<EMAIL_ADDRESS>
###### Abstract

We study the distribution of the matrix product $G_{1}G_{2}\cdots G_{r}$ of
$r$ independent Gaussian matrices of various sizes, where $G_{i}$ is
$d_{i-1}\times d_{i}$, and we denote $p=d_{0}$, $q=d_{r}$, and require
$d_{1}=d_{r-1}$. Here the entries in each $G_{i}$ are standard normal random
variables with mean $0$ and variance $1$. Such products arise in the study of
wireless communication, dynamical systems, and quantum transport, among other
places. We show that, provided each $d_{i}$, $i=1,\ldots,r$, satisfies
$d_{i}\geq Cp\cdot q$, where $C\geq C_{0}$ for a constant $C_{0}>0$ depending
on $r$, then the matrix product $G_{1}G_{2}\cdots G_{r}$ has variation
distance at most $\delta$ to a $p\times q$ matrix $G$ of i.i.d. standard
normal random variables with mean $0$ and variance $\prod_{i=1}^{r-1}d_{i}$.
Here $\delta\rightarrow 0$ as $C\rightarrow\infty$. Moreover, we show a
converse for constant $r$ that if
$d_{i}<C^{\prime}\max\\{p,q\\}^{1/2}\min\\{p,q\\}^{3/2}$ for some $i$, then
this total variation distance is at least $\delta^{\prime}$, for an absolute
constant $\delta^{\prime}>0$ depending on $C^{\prime}$ and $r$. This converse
is best possible when $p=\Theta(q)$.

##  1 Introduction

Random matrices play a central role in many areas of theoretical, applied, and
computational mathematics. One particular application is dimensionality
reduction, whereby one often chooses a rectangular random matrix
$G\in\mathbb{R}^{m\times n}$, $m\ll n$, and computes $G\cdot x$ for a fixed
vector $x\in\mathbb{R}^{n}$. Indeed, this is the setting in compressed sensing
and sparse recovery [12], randomized numerical linear algebra [18, 20, 36],
and sketching algorithms for data streams [25]. Often $G$ is chosen to be a
Gaussian matrix, and in particular, an $m\times n$ matrix with entries that
are i.i.d. normal random variables with mean $0$ and variance $1$, denoted by
$N(0,1)$. Indeed, in compressed sensing, such matrices can be shown to satisfy
the Restricted Isometry Property (RIP) [10], while in randomized numerical
linear algebra, in certain applications such as support vector machines [29]
and non-negative matrix factorization [19], their performance is shown to
often outperform that of other sketching matrices.

Our focus in this paper will be on understanding the product of two or more
Gaussian matrices. Such products arise naturally in different applications.
For example, in the over-constrained ridge regression problem
$\min_{x}\|Ax-b\|_{2}^{2}+\lambda\|x\|_{2}^{2}$, the design matrix
$A\in\mathbb{R}^{n\times d}$, $n\gg d$, is itself often assumed to be Gaussian
(see, e.g., [26]). In this case, the “sketch-and-solve” algorithmic framework
for regression [32] would compute $G\cdot A$ and $G\cdot b$ for an $m\times n$
Gaussian matrix $G$ with $m\approx sd_{\lambda}$, where $sd_{\lambda}$ is the
so-called statistical dimension [2], and solve for the $x$ which minimizes
$\|G\cdot Ax-G\cdot b\|_{2}^{2}+\lambda\|x\|_{2}^{2}$. While computing $G\cdot
A$ is slower than computing the corresponding matrix product for other kinds
of sketching matrices $G$, it often has application-specific [29, 19] as well
as statistical benefits [31]. Notice that $G\cdot A$ is the product of two
independent Gaussian matrices, and in particular, $G$ has a small number of
rows while $A$ has a small number of columns – this is precisely the
rectangular case we will study below. Other applications in randomized
numerical linear algebra where the product of two Gaussian matrices arises is
when one computes the product of a Gaussian sketching matrix and Gaussian
noise in a spiked identity covariance model [37].

The product of two or more Gaussian matrices also arises in diverse fields
such as multiple-input multiple-output (MIMO) wireless communication channels
[24]. Indeed, similar to the above regression problem in which one wants to
reconstruct an underlying vector $x$, in such settings one observes the vector
$y=G_{1}\cdots G_{r}\cdot x+\eta,$ where $x$ is the transmitted signal and
$\eta$ is background noise. This setting corresponds to the situation in which
there are $r$ scattering environments separated by major obstacles, and the
dimensions of the $G_{i}$ correspond to the number of “keyholes” [24]. To
determine the mutual information of this channel, one needs to understand the
singular values of the matrix $G_{1}\cdots G_{r}$. If one can show the
distribution of this product is close to that of a Gaussian distribution in
total variation distance, then one can use the wide range of results known for
the spectrum of a single Gaussian matrix (see, e.g., [35]). Other applications
of products of Gaussian matrices include disordered spin chains [11, 3, 15],
stability of large complex dynamical systems [22, 21], symplectic maps and
Hamiltonian mechanics [11, 4, 28], quantum transport in disordered wires [23,
13], and quantum chromodynamics [27]; we refer the reader to [14, 1] for an
overview.

The main question we ask in this work is:

What is the distribution of the product $G_{1}G_{2}\cdots G_{r}$ of $r$
independent Gaussian matrices of various sizes, where $G_{i}$ is
$d_{i-1}\times d_{i}$?

Our main interest in the question above will be when $G_{1}$ has a small
number $p=d_{0}$ of rows, and $G_{r}$ has a small number $q=d_{r}$ of columns.
Despite the large body of work on random matrix theory (see, e.g., [34] for a
survey), we are not aware of any work which attempts to bound the total
variation distance of the entire distribution of $G_{1}G_{2}\cdots G_{r}$ to a
Gaussian distribution itself.

###  1.1 Our Results

Formally, we consider the problem of distinguishing the product of normalized
Gaussian matrices

$A_{r}=\left(\frac{1}{\sqrt{d_{1}}}G_{1}\right)\left(\frac{1}{\sqrt{d_{2}}}G_{2}\right)\cdots\left(\frac{1}{\sqrt{d_{r-1}}}G_{r-1}\right)\left(\frac{1}{\sqrt{d_{1}}}G_{r}\right)$

from a single normalized Gaussian matrix

$A_{1}=\frac{1}{\sqrt{d_{1}}}G_{1}.$

We show that, when $r$ is a constant, with constant probability we cannot
distinguish the distributions of these two random matrices when $d_{i}\gg p,q$
for all $i$; and, conversely, with constant probability, we can distinguish
these two distributions when the $d_{i}$ are not large enough.

######  Theorem 1 (Main theorem).

Suppose that $d_{i}\geq\max\\{p,q\\}$ for all $i$ and $d_{r-1}=d_{1}$.

  1. (a)

It holds that

$d_{TV}(A_{r},A_{1})\leq C_{1}\sum_{i=1}^{r-1}\sqrt{\frac{pq}{d_{i}}},$

where $d_{TV}(A_{r},A_{1})$ denotes the total variation distance between
$A_{r}$ and $A_{1}$, and $C_{1}>0$ is an absolute constant.

  2. (b)

If $p,q,d_{1},\dots,d_{r}$ further satisfy that

$\sum_{j=1}^{r-1}\frac{1}{d_{j}}\geq\frac{C_{2}^{r}}{\max\\{p,q\\}^{\frac{1}{2}}\min\\{p,q\\}^{\frac{3}{2}}},$

where $C_{2}>0$ is an absolute constant, then $d_{TV}(A_{r},A_{1})\geq 2/3$.

Part (a) states that $d_{TV}(A_{r},A_{1})<2/3$ when $d_{i}\geq
C_{1}^{\prime}pq$ for all $i$ for a constant $C_{1}^{\prime}$ depending on
$r$. The converse in (b) implies that $d_{TV}(A_{r},A_{1})\geq 2/3$ when
$d_{i}\leq C_{2}^{\prime}\max\\{p,q\\}^{1/2}\min\\{p,q\\}^{3/2}$ for some $i$
for a constant $C_{2}^{\prime}$ depending on $r$. When $p=\Theta(q)$ and $r$
is a constant, we obtain a dichotomy (up to a constant factor) for the
conditions on $p,q$ and $d_{i}$.

###  1.2 Our Techniques

#### Upper Bound.

We start by explaining our main insight as to why the distribution of a
product $G_{1}\cdot G_{2}$ of a $p\times d$ matrix $G_{1}$ of i.i.d. $N(0,1)$
random variables and a $d\times q$ matrix $G_{2}$ of i.i.d. $N(0,1)$ random
variables has low variation distance to the distribution of a $p\times q$
matrix $A$ of i.i.d. $N(0,d)$ random variables. One could try to directly
understand the probability density function as was done in the case of Wishart
matrices in [7, 30], which corresponds to the setting when $G_{1}=G_{2}$.
However, there are certain algebraic simplifications in the case of the
Wishart distribution that seem much less tractable when manipulating the
density function of the product of independent Gaussians [9]. Another approach
would be to try to use entropic methods as in [8, 6]. Such arguments try to
reveal entries of the product $G_{1}\cdot G_{2}$ one-by-one, arguing that for
most conditionings of previous entries, the new entry still looks like an
independent Gaussian. However, the entries are clearly not independent – if
$(G_{1}\cdot G_{2})_{i,j}$ has large absolute value, then $(G_{1}\cdot
G_{2})_{i,j^{\prime}}$ is more likely to be large in absolute value, as it
could indicate that the $i$-th row of $G_{1}$ has large norm. One could try to
first condition on the norms of all rows of $G_{1}$ and columns of $G_{2}$,
but additional issues arise when one looks at submatrices: if $(G_{1}\cdot
G_{2})_{i,j},(G_{1}\cdot G_{2})_{i,j^{\prime}},$ and $(G_{1}\cdot
G_{2})_{i^{\prime},j}$ are all large, then it could mean the $i$-th row of
$G_{1}$ and the $i^{\prime}$-th row of $G_{1}$ are correlated with each other,
since they both are correlated with the $j$-th column of $G_{2}$.
Consequently, since $(G_{1}\cdot G_{2})_{i,j^{\prime}}$ is large, it could
make it more likely that $(G_{1}\cdot G_{2})_{i^{\prime},j^{\prime}}$ has
large absolute value. This makes the entropic method difficult to apply in
this context.

Our upper bound instead leverages beautiful work of Jiang [16] and Jiang and
Ma [17] which bounds the total variation distance between the distribution of
an $r\times\ell$ submatrix of a random $d\times d$ orthogonal matrix
(orthonormal rows and columns) and an $r\times\ell$ matrix with i.i.d.
$N(0,1/d)$ entries. Their work shows that if $r\cdot\ell/d\to 0$ as
$d\to\infty$, then the total variation distance between these two matrix
ensembles goes to $0$. It is not immediately clear how to apply such results
in our context. First of all, which submatrix should we be looking at? Note
though, that if $V^{T}$ is a $p\times d$ uniformly random (Haar measure)
matrix with orthonormal rows, and $E$ is a $d\times q$ uniformly random matrix
with orthonormal columns, then by rotational invariance, $V^{T}E$ is
identically distributed to a $p\times q$ submatrix of a $d\times d$ random
orthonormal matrix. Thus, setting $r=p$ and $\ell=q$ in the above results,
they imply that $V^{T}E$ is close in variation distance to a $p\times q$
matrix $H$ with i.i.d. $N(0,1/d)$ entries. Given $G_{1}$ and $G_{2}$, one
could then write them in their singular value decomposition, obtaining
$G_{1}=U\Sigma V^{T}$ and $G_{2}=E\mathrm{T}F^{T}$. Then $V^{T}$ and $E$ are
independent and well-known to be uniformly random $p\times d$ and $d\times q$
orthonormal matrices, respectively. Thus $G_{1}\cdot G_{2}$ is close in total
variation distance to $U\Sigma H\mathrm{T}F^{T}$. However, this does not
immediately help either, as it is not clear what the distribution of this
matrix is. Instead, the “right” way to utilize the results above is to (1)
observe that $G_{1}\cdot G_{2}=U\Sigma V^{T}G_{2}$ is identically distributed
as $U\Sigma X$, where $X$ is a matrix of i.i.d. normal random variables, given
the rotational invariance of the Gaussian distribution. Then (2) $X$ is itself
close to a product $W^{T}Z$ where $W^{T}$ is a random $p\times d$ matrix with
orthonormal rows, and $Z$ is a random $d\times q$ matrix with orthonormal
columns, by the above results. Thus, $G_{1}\cdot G_{2}$ is close to $U\Sigma
W^{T}Z$. Then (3) $U\Sigma W^{T}$ has the same distribution as $G_{1}$, so
$U\Sigma W^{T}Z$ is close to $G_{1}^{\prime}Z$, where $G_{1}^{\prime}$ and
$G_{1}$ are identically distributed, and $G_{1}^{\prime}$ is independent of
$Z$. Finally, (4) $G_{1}^{\prime}Z$ is identically distributed as a matrix
$A_{1}$ of standard normal random variables because $G_{1}^{\prime}$ is
Gaussian and $Z$ has orthonormal columns, by rotational invariance of the
Gaussian distribution.

We hope that this provides a general method for arguments involving Gaussian
matrices - in step (2) we had the quantity $U\Sigma X$, where $X$ was a
Gaussian matrix, and then viewed $X$ as a product of a short-fat random
orthonormal matrix $W^{T}$ and a tall-thin random orthonormal matrix $Z$. Our
proof for the product of more than $2$ matrices recursively uses similar
ideas, and bounds the growth in variation distance as a function of the number
$r$ of matrices involved in the product.

#### Lower Bound.

For our lower bound for constant $r$, we show that the fourth power of the
Schatten 4-norm of a matrix, namely,
$\|X\|_{S_{4}}^{4}=\operatorname{tr}((X^{T}X)^{2})$, can be used to
distinguish a product $A_{r}$ of $r$ Gaussian matrices and a single Gaussian
matrix $A_{1}$. We use Chebyshev’s inequality, for which we need to find the
expectation and variance of $\operatorname{tr}((X^{T}X)^{2})$ for $X=A_{r}$
and $X=A_{1}$.

Let us consider the expectation first. An idea is to calculate the expectation
recursively, that is, for a fixed matrix $M$ and a Gaussian random matrix $G$
we express $\operatorname*{\mathbb{E}}\operatorname{tr}(((MG)^{T}(MG))^{2})$
in terms of $\operatorname*{\mathbb{E}}\operatorname{tr}((M^{T}M)^{2})$. The
real situation turns out to be slightly more complicated. Instead of
expressing $\operatorname*{\mathbb{E}}\operatorname{tr}(((MG)^{T}(MG))^{2})$
in terms of $\operatorname*{\mathbb{E}}\operatorname{tr}((M^{T}M)^{2})$
directly, we decompose
$\operatorname*{\mathbb{E}}\operatorname{tr}(((MG)^{T}(MG))^{2})$ into the sum
of expectations of a few functions in terms of $M$, say,

$\operatorname*{\mathbb{E}}\operatorname{tr}(((MG)^{T}(MG))^{2})=\operatorname*{\mathbb{E}}f_{1}(M)+\operatorname*{\mathbb{E}}f_{2}(M)+\cdots+\operatorname*{\mathbb{E}}f_{s}(M)$

and build up the recurrence relations for
$\operatorname*{\mathbb{E}}f_{1}(MG),\dots,\operatorname*{\mathbb{E}}f_{s}(MG)$
in terms of $\operatorname*{\mathbb{E}}f_{1}(M)$,
$\operatorname*{\mathbb{E}}f_{2}(M)$, …, $\operatorname*{\mathbb{E}}f_{s}(M)$.
It turns out that the recurrence relations are all linear, i.e.,

$\operatorname*{\mathbb{E}}f_{i}(MG)=\sum_{j=1}^{s}a_{ij}\operatorname*{\mathbb{E}}f_{j}(M),\quad
i=1,\dots,s,$ (1)

whence we can solve for $\operatorname*{\mathbb{E}}f_{i}(A_{r})$ and obtaining
the desired expectation
$\operatorname*{\mathbb{E}}\operatorname{tr}((A_{r}^{T}A_{r})^{2})$.

Now we turn to variance. One could try to apply the same idea of finding
recurrence relations for
$\operatorname*{Var}(Q)=\operatorname*{\mathbb{E}}(Q^{2})-(\operatorname*{\mathbb{E}}Q)^{2}$
(where $Q=\operatorname{tr}(((MG)^{T}(MG))^{2})$), but it quickly becomes
intractable for the $\operatorname*{\mathbb{E}}(Q^{2})$ term as it involves
products of eight entries of $M$, which all need to be handled carefully as to
avoid any loose bounds; note, the subtraction of
$(\operatorname*{\mathbb{E}}Q)^{2}$ is critically needed to obtain a small
upper bound on $\operatorname*{Var}(Q)$ and thus loose bounds on
$\operatorname*{\mathbb{E}}(Q^{2})$ would not suffice. For a tractable
calculation, we keep the product of entries of $M$ to $4$th order throughout,
without involving any terms of $8$th order. To do so, we invoke the law of
total variance,

$\operatorname*{Var}_{M,G}(\operatorname{tr}((MG)^{T}(MG))^{2})=\operatorname*{\mathbb{E}}_{M}\left(\left.\operatorname*{Var}_{G}(\operatorname{tr}((G^{T}M^{T}MG)^{2}))\right|M\right)+\operatorname*{Var}_{M}\left(\left.\operatorname*{\mathbb{E}}_{G}\operatorname{tr}((G^{T}M^{T}MG)^{2})\right|M\right).$
(2)

For the first term on the right-hand side, we use Poincaré’s inequality to
upper bound it. Poincaré’s inequality for the Gaussian measure states that for
a differentiable function $f$ on $\mathbb{R}^{m}$,

$\operatorname*{Var}_{g\sim N(0,I_{m})}(f(g))\leq
C\operatorname*{\mathbb{E}}_{g\sim N(0,I_{m})}\|\nabla f(g)\|_{2}^{2}.$

Here we can simply let $f(X)=\operatorname{tr}((MX)^{T}(MX))^{2})$ and
calculate $\operatorname*{\mathbb{E}}\|\nabla f(G)\|_{2}^{2}$. This is
tractable since $\operatorname*{\mathbb{E}}\|\nabla f(G)\|_{2}^{2}$ involves
the products of at most $4$ entries of $M$, and we can use the recursive idea
for the expectation above to express

$\operatorname*{\mathbb{E}}\|\nabla
f(G)\|_{2}^{2}=\sum_{i}a_{ij}\operatorname*{\mathbb{E}}g_{i}(M)$

for a few functions $g_{i}$’s and establish a recurrence relation for each
$g_{i}$.

The second term on the right-hand side of (2) can be dealt with by plugging in
(1), and turns out to depend on a new quantity
$\operatorname*{Var}(\operatorname{tr}^{2}(M^{T}M))$. We again apply the
recursive idea and the law of total variance to

$\operatorname*{Var}_{M,G}(\operatorname{tr}^{2}(G^{T}M^{T}MG))=\operatorname*{\mathbb{E}}_{M}\left(\left.\operatorname*{Var}_{G}(\operatorname{tr}^{2}((G^{T}M^{T}MG))\right|M\right)+\operatorname*{Var}_{M}\left(\left.\operatorname*{\mathbb{E}}_{G}\operatorname{tr}^{2}(G^{T}M^{T}MG)\right|M\right).$

Again, the first term on the right-hand side can be handled by Poincaré’s
inequality and the second-term turns out to depend on
$\operatorname*{Var}(\operatorname{tr}((M^{T}M)^{2}))$, which is crucial. We
have now obtained a double recurrence involving inequalities on
$\operatorname*{Var}(\operatorname{tr}((M^{T}M)^{2}))$ and
$\operatorname*{Var}(\operatorname{tr}^{2}((M^{T}M)^{2}))$, from which we can
solve for an upper bound on
$\operatorname*{Var}(\operatorname{tr}(A_{r}^{T}A_{r})^{2})$. This upper
bound, however, grows exponentially in $r$, which is impossible to improve due
to our use of Poincaré’s inequality.

##  2 Preliminaries

#### Notation.

For a random variable $X$ and a probability distribution $\mathcal{D}$, we use
$X\sim\mathcal{D}$ to denote that $X$ is subject to $\mathcal{D}$. For two
random variables $X$ and $Y$ defined on the same sample space, we write
$X\mathop{\overset{\text{ }d\text{
}}{\leavevmode\resizebox{0.0pt}{0.0pt}{=}}}Y$ if $X$ and $Y$ are identically
distributed.

We use $\mathcal{G}_{m,n}$ to denote the distribution of $m\times n$ Gaussian
random matrices of i.i.d. entries $N(0,1)$ and $\mathcal{O}_{m,n}$ to denote
the uniform distribution (Haar) of an $m\times n$ random matrix with
orthonormal rows. For a distribution $\mathcal{D}$ on a linear space and a
scaling factor $\alpha\in\mathbb{R}$, we use $\alpha\mathcal{D}$ to denote the
distribution of $\alpha X$, where $X\sim\mathcal{D}$.

For two probability measures $\mu$ and $\nu$ on the Borel algebra
$\mathcal{F}$ of $\mathbb{R}^{m}$, the total variation distance between $\mu$
and $\nu$ is defined as

$d_{TV}(\mu,\nu)=\sup_{A\in\mathcal{F}}|\mu(A)-\nu(A)|=\frac{1}{2}\int_{\mathbb{R}^{m}}\left|\frac{d\mu}{d\nu}-1\right|d\nu.$

If $\nu$ is absolutely continuous with respect to $\mu$, one can define the
Kullback-Leibler Divergence between $\mu$ and $\nu$ as

$D_{\mathrm{KL}}(\mu\|\nu)=\int_{\mathbb{R}^{m}}\frac{d\mu}{d\nu}\log_{2}\frac{d\mu}{d\nu}d\nu.$

If $\nu$ is not absolutely continuous with respect to $\mu$, we define
$D_{\mathrm{KL}}(\mu\|\nu)=\infty$.

When $\mu$ and $\nu$ correspond to two random variables $X$ and $Y$,
respectively, we also write $d_{TV}(\mu,\nu)$ and $D_{\mathrm{KL}}(\mu\|\nu)$
as $d_{TV}(X,Y)$ and $D_{\mathrm{KL}}(X\|Y)$, respectively.

The following is the well-known relation between the Kullback-Leibler
divergence and the total variation distance between two probability measures.

######  Lemma 1 (Pinsker’s Inequality [5, Theorem 4.19]).

$d_{TV}(\mu,\nu)\leq\sqrt{\frac{1}{2}D_{\mathrm{KL}}(\mu\|\nu)}$.

The following result, concerning the distance between the submatrix of a
properly scaled Gaussian random matrix and a submatrix of a random orthogonal
matrix, is due to Jiang and Ma [17].

######  Lemma 2 ([17]).

Let $G\sim\mathcal{G}_{d,d}$ and $Z\sim\mathcal{O}_{d,d}$. Suppose that
$p,q\leq d$ and $\hat{G}$ is the top-left $p\times q$ block of $G$ and
$\hat{Z}$ the top-left $p\times q$ block of $Z$. Then

$d_{\mathrm{KL}}\left(\left.\frac{1}{\sqrt{d}}\hat{G}\right\|\hat{Z}\right)\leq
C\frac{pq}{d},$ (3)

where $C>0$ is an absolute constant.

The original paper [17] does not state explicitly the bound in (3) and only
states that the Kullback-Leibler divergence tends to $0$ as $d\to\infty$. A
careful examination of the proof of [17, Theorem 1(i)], by keeping track of
the order of the various $o(1)$ terms, reveals the quantitative bound (3).

#### Useful Inequalities.

We list two useful inequalities below.

######  Lemma 3 (Poincaré’s inequality for Gaussian measure [5, Theorem
3.20]).

Let $X\sim N(0,I_{n})$ be the standard $n$-dimensional Gaussian distribution
and $f:\mathbb{R}^{n}\to\mathbb{R}$ be any continuously differentiable
function. Then

$\operatorname*{Var}(f(X))\leq\operatorname*{\mathbb{E}}\left(\left\|\nabla
f(X)\right\|_{2}^{2}\right).$

######  Lemma 4 (Trace inequality, [33]).

Let $A$ and $B$ be symmetric, positive semidefinite matrices and $k$ be a
positive integer. Then

$\operatorname{tr}((AB)^{k})\leq\min\left\\{\left\|A\right\|_{op}^{k}\operatorname{tr}(B^{k}),\left\|B\right\|_{op}^{k}\operatorname{tr}(A^{k})\right\\}.$

##  3 Upper Bound

Let $r\geq 2$ be an integer. Suppose that $G_{1},\dots,G_{r}$ are independent
Gaussian random matrices, where $G_{i}\sim\mathcal{G}_{d_{i-1},d_{i}}$ and
$d_{0}=p$, $d_{r}=q$ and $d_{r-1}=d_{1}$. Consider the product of normalized
Gaussian matrices

$\displaystyle
A_{r}=\left(\frac{1}{\sqrt{d_{1}}}G_{1}\right)\left(\frac{1}{\sqrt{d_{2}}}G_{2}\right)\cdots\left(\frac{1}{\sqrt{d_{r-1}}}G_{r-1}\right)\left(\frac{1}{\sqrt{d_{1}}}G_{r}\right)$
and a single normalized Gaussian random matrix $\displaystyle
A_{1}=\frac{1}{\sqrt{d_{1}}}G_{1}^{\prime}$

where $G_{1}^{\prime}\sim\mathcal{G}_{p,q}$. In this section, we shall show
that when $p,q\ll d_{i}$ for all $i$, we cannot distinguish $A_{r}$ from
$A_{1}$ with constant probability.

For notational convenience, let $W_{i}=\frac{1}{\sqrt{d_{i}}}G_{i}$ for $i\leq
r$ and $W_{r}=\frac{1}{\sqrt{d_{1}}}G_{r}$. Assume that $pq\leq\beta d_{i}$
for some constant $\beta$ for all $i$. Our question is to find the total
variation distance between the matrix product $W_{1}W_{2}\cdots W_{r}$ and the
product $W_{1}W_{r}$ of two matrices.

######  Lemma 5.

Let $p,q,d,d^{\prime}$ be positive integers satisfying that $pq\leq\beta d$
and $pq\leq\beta d^{\prime}$ for some constant $\beta<1$. Suppose that
$A\in\mathbb{R}^{p\times d}$,
$G\sim\frac{1}{\sqrt{d}}\mathcal{G}_{d,d^{\prime}}$, and
$L\sim\mathcal{O}_{d^{\prime},d}$. Further suppose that $G$ and $L$ are
independent. Let $Z\sim\mathcal{O}_{q,d}$ be independent of $A$, $G$ and $L$.
Then

$d_{TV}(AGL,AZ^{T})\leq C\sqrt{\frac{pq}{d}},$

where $C>0$ is an absolute constant.

###### Proof.

Let $A=U\Sigma V^{T}$ be its singular value decomposition, where $V$ has
dimension $d\times p$. Then

$AGL=U\Sigma(V^{T}GL)\mathop{\overset{\text{ }d\text{
}}{\leavevmode\resizebox{0.0pt}{0.0pt}{=}}}U\Sigma X,$

where $X$ is a $p\times q$ random matrix of i.i.d. $N(0,1/d)$ entries. Suppose
that $\tilde{Z}$ consists of the top $p$ rows of $Z^{T}$. Then

$AZ^{T}=U\Sigma(V^{T}Z^{T})\mathop{\overset{\text{ }d\text{
}}{\leavevmode\resizebox{0.0pt}{0.0pt}{=}}}U\Sigma\tilde{Z}.$

Note that $X$ and $Z$ are independent of $U$ and $\Sigma$. It follows from
Lemma 2 that

$d_{\mathrm{KL}}(AGL\|AZ^{T})=d_{\mathrm{KL}}(U\Sigma
X\|U\Sigma\tilde{Z})=d_{\mathrm{KL}}(X\|\tilde{Z})\leq C\frac{pq}{d},$

where $C>0$ is an absolute constant. The result follows from Pinsker’s
inequality (Lemma 1). ∎

The next theorem follows from the above.

######  Theorem 2.

It holds that

$d_{TV}(W_{1}\cdots W_{r},W_{1}W_{r})\leq
C\sum_{i=1}^{r}\sqrt{\frac{pq}{d_{i}}},$

where $C>0$ is an absolute constant.

###### Proof.

Let $W_{r}=U\Sigma V^{T}$ and $X_{i}\sim\mathcal{O}_{q,d_{i}}$, independent
from each other and from the $W_{i}$’s. Applying the preceding lemma with
$A=W_{1}\cdots W_{r-2}$, $G=W_{r-1}$ and $L=U$, we have

$d_{TV}(W_{1}\cdots W_{r-2}W_{r-1}W_{r},W_{1}\cdots W_{r-2}X_{r-1}^{T}\Sigma
V^{T})\leq C\sqrt{\frac{pq}{d_{r-1}}},$

Next, applying the preceding lemma with $A=W_{1}\cdots W_{r-3}$, $G=W_{r-1}$
and $L=X_{r}$, we have

$d_{TV}(W_{1}\cdots W_{r-2}X_{r}^{T}\Sigma V^{T},W_{1}\cdots
W_{r-3}X_{r-2}^{T}\Sigma V^{T})\leq C\sqrt{\frac{pq}{d_{r-2}}},$

Iterating this procedure, we have in the end that

$d_{TV}(W_{1}W_{2}X_{3}^{T}\Sigma V^{T},W_{1}X_{2}^{T}\Sigma V^{T})\leq
C\sqrt{\frac{pq}{d_{2}}}.$

Since $U$, $\Sigma$ and $V$ are independent and $X_{2}\mathop{\overset{\text{
}d\text{ }}{\leavevmode\resizebox{0.0pt}{0.0pt}{=}}}U^{T}$, it holds that
$X_{2}^{T}\Sigma V^{T}\mathop{\overset{\text{ }d\text{
}}{\leavevmode\resizebox{0.0pt}{0.0pt}{=}}}W_{r}$. Therefore,

$d_{TV}(W_{1}\cdots W_{r},W_{1}W_{r})\leq
C\sum_{i=2}^{r-1}\sqrt{\frac{pq}{d_{i}}}.\qed$

Repeating the same argument for $W_{1}W_{r}$, we obtain the following
corollary immediately.

######  Corollary 1.

It holds that

$d_{TV}(A_{r},A_{1})\leq C\sum_{i=1}^{r-1}\sqrt{\frac{pq}{d_{i}}},$

where $C>0$ is an absolute constant.

##  4 Lower Bound

Suppose that $r$ is a constant. We shall show that one can distinguish the
product of $r$ Gaussian random matrices

$A_{r}=\left(\frac{1}{\sqrt{d_{1}}}G_{1}\right)\left(\frac{1}{\sqrt{d_{2}}}G_{2}\right)\cdots\left(\frac{1}{\sqrt{d_{r-1}}}G_{r-1}\right)\left(\frac{1}{\sqrt{d_{1}}}G_{r}\right),$

from one Gaussian random matrix

$A_{1}=\frac{1}{\sqrt{d_{1}}}G_{1}^{\prime}$

when the intermediate dimensions $d_{1},\dots,d_{r-1}$ are not large enough.
Considering $h(X)=\operatorname{tr}((X^{T}X)^{2})$, it suffices to show that
one can distinguish $h(A_{r})$ and $h(A_{1})$ with a constant probability for
constant $r$. By Chebyshev’s inequality, it suffices to show that

$\max\left\\{\sqrt{\operatorname*{Var}(h(A_{1}))},\sqrt{\operatorname*{Var}(h(A_{r}))}\right\\}\leq
c(\operatorname*{\mathbb{E}}h(A_{r})-\operatorname*{\mathbb{E}}h(A_{1}))$

for a small constant $c$. We calculate that:

######  Lemma 6.

Suppose that $r$ is a constant, $d_{i}\geq\max\\{p,q\\}$ for all
$i=1,\dots,r$. When $p,q,d_{1},\dots,d_{r}\to\infty$,

$\operatorname*{\mathbb{E}}h(A_{r})=\frac{pq(p+q+1)}{d_{r}^{2}}+(1+o(1))\frac{pq(p-1)(q-1)}{d_{r}^{2}}\sum_{j=1}^{r-1}\frac{1}{d_{j}}.$

######  Lemma 7.

Suppose that $r$ is a constant, $d_{i}\geq\max\\{p,q\\}$ for all
$i=1,\dots,r$. There exists an absolute constant $C$ such that, when
$p,q,d_{1},\dots,d_{r}$ are sufficently large,

$\operatorname*{Var}(h(A_{r}))\leq\frac{C^{r}(p^{3}q+pq^{3})}{d_{1}^{4}}.$

We conclude with the following theorem, which can be seen as a tight converse
to Corollary 1 up to a constant factor on the conditions for
$p,q,d_{1},\dots,d_{r}$.

######  Theorem 3.

Suppose that $r$ is a constant and $d_{i}\geq\max\\{p,q\\}$ for all
$i=1,\dots,r$. Further suppose that $d_{1}=d_{r}$. When
$p,q,d_{1},\dots,d_{r}$ are sufficiently large and satisfy that

$\sum_{j=1}^{r-1}\frac{1}{d_{j}}\geq\frac{C^{r}}{\max\\{p,q\\}^{\frac{1}{2}}\min\\{p,q\\}^{\frac{3}{2}}},$

where $C>0$ is some absolute constant, with probability at least $2/3$, one
can distinguish $A_{r}$ from $A_{1}$.

###  4.1 Calculation of the Mean

Suppose that $A$ is a $p\times q$ random matrix, and is rotationally invariant
under left- and right-multiplication by orthogonal matrices. We define

$\displaystyle S_{1}(p,q)$
$\displaystyle=\operatorname*{\mathbb{E}}A_{11}^{4}\quad\text{(diagonal)}$
$\displaystyle S_{2}(p,q)$
$\displaystyle=\operatorname*{\mathbb{E}}A_{21}^{4}\quad\text{(off-diagonal)}$
$\displaystyle S_{3}(p,q)$
$\displaystyle=\operatorname*{\mathbb{E}}A_{i1}^{2}A_{j1}^{2}\quad(i\neq
j)\quad\text{(same column)}$ $\displaystyle S_{4}(p,q)$
$\displaystyle=\operatorname*{\mathbb{E}}A_{1i}^{2}A_{1j}^{2}\quad(i\neq
j)\quad\text{(same row)}$ $\displaystyle S_{5}(p,q)$
$\displaystyle=\operatorname*{\mathbb{E}}A_{1i}^{2}A_{2j}^{2}\quad(i\neq j)$
$\displaystyle S_{6}(p,q)$
$\displaystyle=\operatorname*{\mathbb{E}}A_{ik}A_{il}A_{jk}A_{jl}\quad(i\neq
j,k\neq l)\quad\text{(rectangle)}$

Since $A$ is left- and right-invariant under rotations, these quantities are
well-defined. Then

$\displaystyle\begin{aligned}
\operatorname*{\mathbb{E}}\operatorname{tr}((A^{T}A)^{2})=\operatorname*{\mathbb{E}}\sum_{1\leq
i,j\leq
q}(A^{T}A)_{ij}^{2}&=\sum_{i=1}^{q}\operatorname*{\mathbb{E}}(A^{T}A)_{ii}^{2}+\sum_{1\leq
i,j\leq q,i\neq j}\operatorname*{\mathbb{E}}(A^{T}A)_{ij}^{2}\\\
&=q\operatorname*{\mathbb{E}}(A^{T}A)_{11}^{2}+q(q-1)\operatorname*{\mathbb{E}}(A^{T}A)_{12}^{2}\end{aligned}$
and $\displaystyle\begin{aligned}
\operatorname*{\mathbb{E}}(A^{T}A)_{11}^{2}=\operatorname*{\mathbb{E}}\bigg{(}\sum_{i=1}^{p}A_{i1}^{2}\bigg{)}^{2}&=\sum_{i=1}^{p}\operatorname*{\mathbb{E}}A_{i1}^{4}+\sum_{1\leq
i,j\leq p,i\neq j}\operatorname*{\mathbb{E}}A_{i1}^{2}A_{j1}^{2}\\\
&=\operatorname*{\mathbb{E}}A_{11}^{4}+(p-1)\operatorname*{\mathbb{E}}A_{21}^{4}+p(p-1)\operatorname*{\mathbb{E}}A_{11}^{2}A_{21}^{2}\\\
&=:S_{1}(p,q)+(p-1)S_{2}(p,q)+p(p-1)S_{3}(p,q)\end{aligned}$
$\displaystyle\begin{aligned}
\operatorname*{\mathbb{E}}(A^{T}A)_{12}^{2}=\operatorname*{\mathbb{E}}\bigg{(}\sum_{i=1}^{p}A_{i1}A_{i2}\bigg{)}^{2}&=\sum_{i=1}^{p}\operatorname*{\mathbb{E}}A_{i1}^{2}A_{i2}^{2}+\sum_{1\leq
i,j\leq p,i\neq j}\operatorname*{\mathbb{E}}A_{i1}A_{i2}A_{j1}A_{j2}\\\
&=pS_{4}(p,q)+p(p-1)S_{6}(p,q).\end{aligned}$

When $S_{1}(p,q)=S_{2}(p,q)$, we have

$\displaystyle\operatorname*{\mathbb{E}}\operatorname{tr}((A^{T}A)^{2})$
$\displaystyle=q(pS_{1}(p,q)+p(p-1)S_{3}(p,q))+q(q-1)(pS_{4}(p,q)+p(p-1)S_{6}(p,q))$
$\displaystyle=pqS_{1}(p,q)\\!+\\!pq(p\\!-\\!1)S_{3}(p,q)\\!+\\!pq(q\\!-\\!1)S_{4}(p,q)\\!+\\!p(p\\!-\\!1)q(q\\!-\\!1)S_{6}(p,q).$

When $A=G$, we have

$\displaystyle S_{1}(p,q)=S_{2}(p,q)=3,\quad
S_{3}(p,q)=S_{4}(p,q)=S_{5}(p,q)=1,\quad S_{6}(p,q)=0$ and so
$\displaystyle\operatorname*{\mathbb{E}}\operatorname{tr}((A^{T}A)^{2})=3pq+pq(p-1)+pq(q-1)=pq(p+q+1).$

Next, consider $A=BG$, where $B$ is a $p\times d$ random matrix and $G$ a
$d\times q$ random matrix of i.i.d. $N(0,1)$ entries. The following
proposition is easy to verify, and its proof is postponed to Appendix
LABEL:sec:proof_of_simple_prop.

######  Proposition 1.

It holds that
$\operatorname*{\mathbb{E}}A_{21}^{4}=\operatorname*{\mathbb{E}}A_{11}^{4}$.

###### Proof.

We have

$\displaystyle\begin{aligned}
\operatorname*{\mathbb{E}}A_{11}^{4}=\operatorname*{\mathbb{E}}\bigg{(}\sum_{i}B_{1i}G_{i1}\bigg{)}^{4}&=\sum_{i,j,k,l}\operatorname*{\mathbb{E}}B_{1i}B_{1j}B_{1k}B_{1l}\operatorname*{\mathbb{E}}G_{i1}G_{j1}G_{k1}G_{l1}\\\
&=3\sum_{i}\operatorname*{\mathbb{E}}B_{1i}^{4}+3\sum_{i\neq
j}\operatorname*{\mathbb{E}}B_{1i}^{2}B_{1j}^{2}\\\ \end{aligned}$ and
$\displaystyle\begin{aligned}
\operatorname*{\mathbb{E}}A_{21}^{4}=\operatorname*{\mathbb{E}}\bigg{(}\sum_{i}B_{2i}G_{i1}\bigg{)}^{4}&=\sum_{i,j,k,l}\operatorname*{\mathbb{E}}B_{2i}B_{2j}B_{2k}B_{2l}\operatorname*{\mathbb{E}}G_{i1}G_{j1}G_{k1}G_{l1}\\\
&=3\sum_{i}\operatorname*{\mathbb{E}}B_{2i}^{4}+3\sum_{i\neq
j}\operatorname*{\mathbb{E}}B_{2i}^{2}B_{2j}^{2}=\operatorname*{\mathbb{E}}A_{11}^{4}.\qed\end{aligned}$

Suppose that the associated functions of $B$ are named
$T_{1},T_{2},T_{3},T_{4},T_{6},T_{5}$. Then we can calculate that (detailed
calculations can be found in Appendix A)

$\displaystyle S_{1}(p,q)$ $\displaystyle=3dT_{1}(p,d)+3d(d-1)T_{4}(p,d)$
$\displaystyle S_{3}(p,q)$
$\displaystyle=3dT_{3}(p,d)+d(d-1)T_{5}(p,d)+2d(d-1)T_{6}(p,d)$ $\displaystyle
S_{4}(p,q)$ $\displaystyle=dT_{1}(p,d)+d(d-1)T_{4}(p,d)$ $\displaystyle
S_{5}(p,q)$ $\displaystyle=dT_{3}(p,d)+d(d-1)T_{5}(p,d)$ $\displaystyle
S_{6}(p,q)$ $\displaystyle=dT_{3}(p,d)+d(d-1)T_{6}(p,d)$

It is clear that $S_{1},S_{3},S_{4},S_{5},S_{6}$ depend only on $d$ (not on
$p$ and $q$) if $T_{1},T_{3},T_{4},T_{5},T_{6}$ do so. Furthermore, if
$T_{1}=3T_{4}$ then we have $S_{1}=3S_{4}$ and thus $S_{4}=d(d+2)T_{4}$. If
$T_{3}=2T_{6}+T_{5}$ then $S_{3}=d(d+2)T_{3}$ and $S_{3}=2S_{6}+S_{5}$. Hence,
if $T_{3}=T_{4}$ then $S_{3}=S_{4}$. We can verify that all these conditions
are satisfied with one Gaussian matrix and we can iterate it to obtain these
quantities for the product of $r$ Gaussian matrices with intermediate
dimensions $d_{1},d_{2},\dots,d_{r-1}$. We have that

$S_{3}=S_{4}=\prod_{i=1}^{r-1}d_{i}(d_{i}+2),\quad S_{1}=3S_{4},\quad
S_{6}=\sum_{j=1}^{r-1}\left(\prod_{i=1}^{j-1}d_{i}(d_{i}+2)\right)d_{j}\left(\prod_{i=j+1}^{r-1}d_{i}(d_{i}-1)\right).$

Therefore, normalizing the $i$-th matrix by $1/\sqrt{d_{i}}$, that is,

$A=\left(\frac{1}{\sqrt{d_{1}}}G_{1}\right)\left(\frac{1}{\sqrt{d_{2}}}G_{2}\right)\cdots\left(\frac{1}{\sqrt{d_{r-1}}}G_{r-1}\right)\left(\frac{1}{\sqrt{d_{1}}}G_{r}\right),$

we have for constant $r$ that

$\displaystyle\operatorname*{\mathbb{E}}\operatorname{tr}((A^{T}A)^{2})$
$\displaystyle=\frac{1}{d_{1}^{2}d_{2}^{2}\cdots
d_{r-1}^{2}d_{1}^{2}}\left(pq(p+q+1)S_{3}+pq(p-1)(q-1)S_{6}\right)$ (4)
$\displaystyle\approx\frac{pq(p+q+1)}{d_{r}^{2}}+\frac{pq(p-1)(q-1)}{d_{r}^{2}}\sum_{j=1}^{r-1}\frac{1}{d_{j}}.$

###  4.2 Calculation of the Variance

Let $M\in\mathbb{R}^{p\times p}$ be a random symmetric matrix, and let
$G\in\mathbb{R}^{p\times q}$ be a random matrix of i.i.d. $N(0,1)$ entries. We
want to find the variance of $\operatorname{tr}((G^{T}MG)^{2})$. The detailed
calculations of some steps can be found in Appendix B.

Our starting point is the law of total variance, which states that

$\operatorname*{Var}(\operatorname{tr}((G^{T}MG)^{2}))=\operatorname*{\mathbb{E}}_{M}\left(\left.\operatorname*{Var}_{G}(\operatorname{tr}((G^{T}MG)^{2}))\right|M\right)+\operatorname*{Var}_{M}\left(\left.\operatorname*{\mathbb{E}}_{G}\operatorname{tr}((G^{T}MG)^{2})\right|M\right)$
(5)

#### Step 1a.

We shall handle each term separately. Consider the first term, which we shall
bound using the Poincaré inequality for Gaussian measures. Define
$f(X)=\operatorname{tr}((X^{T}MX)^{2})$, where $X\in\mathbb{R}^{p\times q}$.
We shall calculate $\nabla f$.

$f(X)=\left\|X^{T}MX\right\|_{F}^{2}=\sum_{1\leq i,j\leq
q}(X^{T}MX)_{ij}^{2}=\sum_{1\leq i,j\leq q}\bigg{(}\sum_{1\leq k,l\leq
p}M_{kl}X_{ki}X_{lj}\bigg{)}^{2}.$

Then

$\frac{\partial f}{\partial X_{rs}}=\sum_{1\leq i,j\leq q}2\bigg{(}\sum_{1\leq
u,v\leq p}M_{uv}X_{ui}X_{vj}\bigg{)}\bigg{(}\sum_{1\leq k,l\leq
p}\frac{\partial}{\partial X_{rs}}(M_{kl}X_{ki}X_{lj})\bigg{)}.$

Note that

$\frac{\partial}{\partial
X_{rs}}(M_{kl}X_{ki}X_{lj})=\begin{cases}M_{kl}X_{lj},&(k,i)=(r,s)\text{ and
}(l,j)\neq(r,s)\\\ M_{kl}X_{ki},&(k,i)\neq(r,s)\text{ and }(l,j)=(r,s)\\\
2M_{rr}X_{rs},&(k,i)=(r,s)\text{ and }(l,j)=(r,s)\\\
0,&\text{otherwise}.\end{cases}$

we have that

$\displaystyle\frac{\partial f}{\partial X_{rs}}$
$\displaystyle=4\left(\sum_{1\leq u,v\leq
p}M_{uv}X_{us}X_{vs}\right)M_{rr}X_{rs}+2\sum_{(l,j)\neq(r,s)}\left(\sum_{1\leq
u,v\leq p}M_{uv}X_{us}X_{vj}\right)M_{rl}X_{lj}$
$\displaystyle\qquad+2\sum_{(k,i)\neq(r,s)}\left(\sum_{1\leq u,v\leq
p}M_{uv}X_{ui}X_{vs}\right)M_{kr}X_{ki}$
$\displaystyle=4\left[\left(\sum_{1\leq u,v\leq
p}M_{uv}X_{us}X_{vs}\right)M_{rr}X_{rs}+\sum_{(l,j)\neq(r,s)}\left(\sum_{u,v}M_{uv}X_{us}X_{vj}\right)M_{rl}X_{lj}\right]$
$\displaystyle=4\sum_{l,j}\left(\sum_{u,v}M_{uv}X_{us}X_{vj}\right)M_{rl}X_{lj}.$

Next we calculate $\operatorname*{\mathbb{E}}(\partial f/\partial X_{rs})^{2}$
when $X$ is i.i.d. Gaussian.

$\left(\frac{1}{4}\frac{\partial f}{\partial
X_{rs}}\right)^{2}=\sum_{\begin{subarray}{c}l,j\\\
l^{\prime},j^{\prime}\end{subarray}}\sum_{\begin{subarray}{c}u,v\\\
u^{\prime},v^{\prime}\end{subarray}}M_{uv}M_{u^{\prime}v^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{us}X_{u^{\prime}s}X_{vj}X_{lj}X_{v^{\prime}j^{\prime}}X_{l^{\prime}j^{\prime}}$

We discuss different cases of $j,j^{\prime},s$.

When $j\neq j^{\prime}\neq s$, it must hold that $u=u^{\prime}$, $v=l$ and
$v^{\prime}=l^{\prime}$ for a possible nonzero contribution, and the total
contribution in this case is at most $q(q-1)B^{(1)}_{r,s}$, where

$B^{(1)}_{r,s}=\sum_{1\leq l,l^{\prime}\leq
p}\sum_{u}M_{ul}M_{ul^{\prime}}M_{rl}M_{rl^{\prime}}=\sum_{u}\langle
M_{u,\cdot},M_{r,\cdot}\rangle^{2}.$

When $j=j^{\prime}\neq s$, it must hold that $u=u^{\prime}$ for a possible
nonzero contribution, and the total contribution in this case is at most
$(q-1)B^{(2)}_{r,s}$, where

$\displaystyle B^{(2)}_{r,s}$
$\displaystyle=\sum_{l,l^{\prime}}\sum_{u,v,v^{\prime}}M_{uv}M_{uv^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{us}^{2}X_{vj}X_{lj}X_{v^{\prime}j}X_{l^{\prime}j}$
$\displaystyle=\left\|M\right\|_{F}^{2}\left\|M_{r,\cdot}\right\|_{2}^{2}+2\sum_{u}\langle
M_{u,\cdot},M_{r,\cdot}\rangle^{2}.$

When $j=s\neq j^{\prime}$, it must hold that $v^{\prime}=l^{\prime}$ for
possible nonzero contribution, and the total contribution in this case is at
most $(q-1)B^{(3)}_{r,s}$, where

$\displaystyle B^{(3)}_{r,s}$ $\displaystyle=\sum_{j^{\prime}\neq
s}\left[\sum_{l,l^{\prime}}\sum_{u,v}M_{uv}M_{u^{\prime}l^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{us}X_{u^{\prime}s}X_{vs}X_{ls}X_{l^{\prime}j^{\prime}}^{2}\right]$
$\displaystyle=\sum_{l,l^{\prime}}(2\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle+\operatorname{tr}(M)M_{ll^{\prime}})M_{rl}M_{rl^{\prime}}.$

When $j=j^{\prime}=s$, the nonzero contribution is

$B^{(4)}_{r,s}=\sum_{l,l^{\prime}}\sum_{\begin{subarray}{c}u,v\\\
u^{\prime},v^{\prime}\end{subarray}}M_{uv}M_{u^{\prime}v^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{us}X_{u^{\prime}s}X_{vs}X_{ls}X_{v^{\prime}s}X_{l^{\prime}s}.$

Since $u,u^{\prime},v,v^{\prime},l,l^{\prime}$ needs to be paired, the only
case which is not covered by $B^{(1)}_{rs},B^{(3)}_{rs}$ and $B^{(3)}_{rs}$ is
when $u=v$, $u^{\prime}=v^{\prime}$ and $l=l^{\prime}$, in which case the
contribution is at most

$\sum_{l}\sum_{u,u^{\prime}}M_{uu}M_{u^{\prime}u^{\prime}}M_{rl}^{2}\operatorname*{\mathbb{E}}X_{us}^{2}X_{u^{\prime}s}^{2}X_{ls}^{2}\lesssim\operatorname{tr}^{2}(M)\left\|M_{r,\cdot}\right\|_{2}^{2}.$

Hence

$B^{(4)}_{r,s}\lesssim
B^{(1)}_{r,s}+B^{(2)}_{r,s}+B^{(3)}_{r,s}+\operatorname{tr}^{2}(M)\left\|M_{r,\cdot}\right\|_{2}^{2}.$

It follows that

$\displaystyle\sum_{r,s}B^{(1)}_{r,s}=q\sum_{u,r}\langle
M_{u,\cdot},M_{r,\cdot}\rangle^{2}=q\operatorname{tr}(M^{4})$
$\displaystyle\sum_{r,s}B^{(2)}_{r,s}=q\sum_{r}\left\|M\right\|_{F}^{2}\left\|M_{r,\cdot}\right\|_{2}^{2}+2q\sum_{u,r}\langle
M_{u,\cdot},M_{r,\cdot}\rangle^{2}=q\left\|M\right\|_{F}^{4}+2q\operatorname{tr}(M^{4})$
$\displaystyle\begin{aligned}
\sum_{r,s}B^{(3)}_{r,s}&=\sum_{r,s}\sum_{l,l^{\prime}}(2\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle+\operatorname{tr}(M)M_{l^{\prime}l^{\prime}})M_{rl}M_{rl^{\prime}}\\\
&\leq
2q\operatorname{tr}(M^{4})+q\operatorname{tr}(M)\left\|M\right\|_{F}\sqrt{\operatorname{tr}(M^{4})}\end{aligned}$

Note that
$\operatorname{tr}(M^{4})\leq\operatorname{tr}^{2}(M^{2})=\left\|M\right\|_{F}^{4}$.
Hence

$\displaystyle\frac{1}{16}\operatorname*{\mathbb{E}}\left\|\nabla
f\right\|_{2}^{2}$
$\displaystyle\leq\sum_{r,s}((q-1)(q-2)B_{rs}^{(1)}+(q-1)B_{rs}^{(2)}+(q-1)B_{rs}^{(3)}+B_{rs}^{(4)})$
$\displaystyle\lesssim\sum_{r,s}(q^{2}B_{rs}^{(1)}+qB_{rs}^{(2)}+qB_{rs}^{(3)}+\operatorname{tr}^{2}(M)\left\|M_{r,\cdot}\right\|_{2}^{2})$
$\displaystyle\lesssim
q^{3}\operatorname{tr}(M^{4})+q^{2}\left\|M\right\|_{F}^{4}+q^{2}\operatorname{tr}(M)\left\|M\right\|_{F}\sqrt{\operatorname{tr}(M^{4})}+q\operatorname{tr}^{2}(M)\left\|M\right\|_{F}^{2}.$

By the Gaussian Poincaré inequality,

$\displaystyle\quad\
\operatorname*{Var}_{G}(\operatorname{tr}((G^{T}MG)^{2})|M)$ (6)
$\displaystyle\lesssim\operatorname*{\mathbb{E}}\left\|\nabla
f\right\|_{2}^{2}$ $\displaystyle\lesssim
q^{3}\operatorname{tr}(M^{4})+q^{2}\left\|M\right\|_{F}^{4}+q^{2}\operatorname{tr}(M)\left\|M\right\|_{F}\sqrt{\operatorname{tr}(M^{4})}+q\operatorname{tr}^{2}(M)\left\|M\right\|_{F}^{2}.$

For the terms on the right-hand side, we calculate that (using the trace
inequality (Lemma 4))

$\displaystyle\begin{aligned}
\operatorname*{\mathbb{E}}\operatorname{tr}((G^{T}MG)^{4})=\operatorname*{\mathbb{E}}\operatorname{tr}((MGG^{T})^{4})\leq\operatorname*{\mathbb{E}}\left\|GG^{T}\right\|_{op}^{4}\operatorname{tr}(M^{4})&=\operatorname*{\mathbb{E}}\left\|G\right\|_{op}^{8}\operatorname{tr}(M^{4})\\\
&\lesssim\max\\{p,q\\}^{4}\operatorname{tr}(M^{4}),\end{aligned}$
$\displaystyle\operatorname*{\mathbb{E}}\left\|G^{T}MG\right\|_{F}^{4}\leq\operatorname*{\mathbb{E}}\left\|G\right\|_{op}^{8}\left\|M\right\|_{F}^{4}\lesssim\max\\{p,q\\}^{4}\left\|M\right\|_{F}^{4},$
$\displaystyle\operatorname*{\mathbb{E}}\operatorname{tr}^{2}(G^{T}MG)\left\|G^{T}MG\right\|_{F}^{2}\leq\operatorname*{\mathbb{E}}\left\|G\right\|_{op}^{8}\operatorname{tr}^{2}(M)\left\|M\right\|_{F}^{2}\lesssim\max\\{p,q\\}^{4}\operatorname{tr}^{2}(M)\left\|M\right\|_{F}^{2}$

and

$\displaystyle\quad\
\operatorname*{\mathbb{E}}\operatorname{tr}(G^{T}MG)\left\|G^{T}MG\right\|_{F}\sqrt{\operatorname{tr}((G^{T}MG)^{4})}$
$\displaystyle\leq\operatorname*{\mathbb{E}}\left\|G\right\|_{op}^{2}\operatorname{tr}(G)\cdot\left\|G\right\|_{op}^{2}\left\|M\right\|_{F}^{2}\cdot\sqrt{\left\|G\right\|_{op}^{8}\operatorname{tr}(M^{4})}$
$\displaystyle=\operatorname*{\mathbb{E}}\left\|G\right\|_{op}^{8}\operatorname{tr}(M)\left\|M\right\|_{F}\sqrt{\operatorname{tr}(M^{4})}$
$\displaystyle\lesssim\max\\{p,q\\}^{4}\operatorname{tr}(M)\left\|M\right\|_{F}\sqrt{\operatorname{tr}(M^{4})}.$

This implies that each term on the right-hand of (6) grows geometrically.

#### Step 1b.

Next we deal with the second term in (5). We have

$\displaystyle\operatorname*{\mathbb{E}}_{G}\operatorname{tr}\left((G^{T}MG)^{2}\right)=\sum_{i,j}\operatorname*{\mathbb{E}}_{G}(G^{T}MG)_{ij}^{2}$
$\displaystyle=\sum_{i,j}\operatorname*{\mathbb{E}}_{G}\bigg{(}\sum_{k,l}M_{kl}G_{ki}G_{lj}\bigg{)}^{2}$
$\displaystyle=\sum_{i,j}\sum_{k,l,k^{\prime},l^{\prime}}M_{kl}M_{k^{\prime}l^{\prime}}\operatorname*{\mathbb{E}}_{G}G_{ki}G_{lj}G_{k^{\prime}i}G_{l^{\prime}j}.$

When $i\neq j$, for non-zero contribution, it must hold that $k=l$ and
$k^{\prime}=l^{\prime}$ and thus the nonzero contribution is

$\sum_{i\neq j}\sum_{k,l}M_{kl}^{2}=q(q-1)\left\|M\right\|_{F}^{2}.$

When $i=j$, the contribution is

$\sum_{i}\sum_{\begin{subarray}{c}k,l,k^{\prime},l^{\prime}\end{subarray}}M_{kl}M_{k^{\prime}l^{\prime}}\operatorname*{\mathbb{E}}G_{ki}G_{li}G_{k^{\prime}i}G_{l^{\prime}i}=2q\left\|M\right\|_{F}^{2}+q\operatorname{tr}^{2}(M).$
(7)

Hence

$\operatorname*{\mathbb{E}}_{G}\operatorname{tr}\left((G^{T}MG)^{2}\right)=q(q+1)\left\|M\right\|_{F}^{2}+q\operatorname{tr}^{2}(M)$

and when $M$ is random,

$\displaystyle\quad\,\operatorname*{Var}\left(\left.\operatorname*{\mathbb{E}}\operatorname{tr}((G^{T}MG)^{2})\right|M\right)$
(8)
$\displaystyle=\operatorname*{Var}\left(q(q+1)\left\|M\right\|_{F}^{2}+q\operatorname{tr}^{2}(M)\right)$
$\displaystyle\leq
q^{2}(q+1)^{2}\operatorname*{Var}(\left\|M\right\|_{F}^{2})+q^{2}\operatorname*{Var}(\operatorname{tr}^{2}(M))+2q^{2}(q+1)\sqrt{\operatorname*{Var}(\left\|M\right\|_{F}^{2})\operatorname*{Var}(\operatorname{tr}^{2}(M))}.$

#### Step 2a.

Note that the $\operatorname*{Var}(\operatorname{tr}^{2}(M))$ term on the
right-hand side of (8). To bound this term, we examine the variance of $g(G)$,
where $g(X)=\operatorname{tr}^{2}(X^{T}MX)$. We shall again calculate $\nabla
g$. Note that

$\frac{\partial g}{\partial
X_{rs}}=2\operatorname{tr}(X^{T}MX)\sum_{i}\sum_{k,l}M_{kl}\frac{\partial}{\partial
X_{rs}}X_{ki}X_{li}$

and

$\frac{\partial}{\partial
X_{rs}}(X_{ki}X_{li})=\begin{cases}X_{li},&(k,i)=(r,s)\text{ and
}(l,i)\neq(r,s)\\\ X_{ki},&(k,i)\neq(r,s)\text{ and }(l,i)=(r,s)\\\
2X_{rs},&(k,i)=(r,s)\text{ and }(l,i)=(r,s)\\\
0,&\text{otherwise}.\end{cases}$

We have

$\frac{\partial g}{\partial
X_{rs}}=4\operatorname{tr}(X^{T}MX)\sum_{l}M_{rl}X_{ls}=4\sum_{\begin{subarray}{c}1\leq
j\leq q\\\ 1\leq l,u,v\leq p\end{subarray}}M_{uv}M_{rl}X_{ls}X_{uj}X_{vj}$

Next we calculate $\operatorname*{\mathbb{E}}(\partial g/\partial X_{rs})^{2}$
when $X$ is i.i.d. Gaussian.

$\left(\frac{1}{4}\frac{\partial g}{\partial
X_{rs}}\right)^{2}=\sum_{\begin{subarray}{c}j,l,u,v\\\
j^{\prime},l^{\prime},u^{\prime},v^{\prime}\end{subarray}}M_{uv}M_{u^{\prime}v^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{ls}X_{l^{\prime}s}X_{uj}X_{vj}X_{u^{\prime}j^{\prime}}X_{v^{\prime}j^{\prime}}$

In order for the expectation in the summand to be non-zero, we must have one
of the following cases: (1) $s\neq j\neq j^{\prime}$, (2) $s=j\neq
j^{\prime}$, (3) $s=j^{\prime}\neq j$, (4) $s\neq j=j^{\prime}$, (5)
$s=j=j^{\prime}$. We calculate the contribution in each case below.

Case 1: it must hold that $l=l^{\prime}$, $u=v$ and $u^{\prime}=v^{\prime}$.
The contribution is $(q-1)(q-2)B^{(1)}_{rs}$, where

$B^{(1)}_{rs}=\sum_{l,u,u^{\prime}}M_{uu}M_{u^{\prime}u^{\prime}}M_{rl}^{2}=\operatorname{tr}^{2}(M)\left\|M_{r,\cdot}\right\|_{2}^{2}.$

Case 2: it must hold that $u^{\prime}=v^{\prime}$. The contribution is
$(q-1)B^{(2)}_{rs}$, where

$\displaystyle B^{(2)}_{rs}$
$\displaystyle=\sum_{l,l^{\prime},u,u^{\prime},v}M_{uv}M_{u^{\prime}u^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{ls}X_{l^{\prime}s}X_{us}X_{vs}X_{u^{\prime}j^{\prime}}^{2}$
$\displaystyle=\operatorname{tr}(M)\bigg{(}\operatorname{tr}(M)\left\|M_{r,\cdot}\right\|_{2}^{2}+2\sum_{l,l^{\prime}}M_{ll^{\prime}}M_{rl}M_{rl^{\prime}}\bigg{)}$

Case 3: this gives the same bound as Case 2.

Case 4: it must hold that $l=l^{\prime}$. The contribution is
$(q-1)B^{(4)}_{rs}$, where

$B^{(4)}_{rs}=\sum_{l,u,u^{\prime},v,v^{\prime}}M_{uv}M_{u^{\prime}v^{\prime}}M_{rl}^{2}\operatorname*{\mathbb{E}}X_{uj}X_{vj}X_{u^{\prime}j}X_{v^{\prime}j}=3\left\|M_{r,\cdot}\right\|_{2}^{2}\left\|M\right\|_{F}^{2}$

Case 5: the contribution is $B^{(5)}_{rs}$, where

$B^{(5)}_{rs}=\sum_{\begin{subarray}{c}l,u,v\\\
l^{\prime},u^{\prime},v^{\prime}\end{subarray}}M_{uv}M_{u^{\prime}v^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{ls}X_{us}X_{vs}X_{l^{\prime}s}X_{u^{\prime}s}X_{v^{\prime}s}.$

The only uncovered case is $l=u^{\prime}$, $l^{\prime}=v$, $u=v^{\prime}$ and
its symmetries. In such a case the contribution is at most

$C\sum_{l,u,v}M_{uv}M_{lu}M_{rl}M_{rv}=C\sum_{u}\langle
M_{r,\cdot},M_{u,\cdot}\rangle^{2}.$

Note that

$\displaystyle\sum_{r,s}B_{rs}^{(1)}$
$\displaystyle=q\operatorname{tr}^{2}(M)\left\|M\right\|_{F}^{2},$
$\displaystyle\sum_{r,s}B_{rs}^{(2)}$
$\displaystyle=q\operatorname{tr}^{2}(M)\left\|M\right\|_{F}^{2}+2q\operatorname{tr}(M)\sum_{l,l^{\prime}}M_{ll^{\prime}}\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle$ $\displaystyle\leq
q\operatorname{tr}^{2}(M)\left\|M\right\|_{F}^{2}+2q\operatorname{tr}(M)\left\|M\right\|_{F}\sqrt{\operatorname{tr}(M^{4})},$
$\displaystyle\sum_{r,s}B_{rs}^{(4)}$
$\displaystyle=q\left\|M\right\|_{F}^{4},$
$\displaystyle\sum_{r,s}B_{rs}^{(5)}$
$\displaystyle\lesssim\sum_{r,s}B_{rs}^{(1)}+\sum_{r,s}B_{rs}^{(2)}+\operatorname{tr}(M^{4}).$

Therefore,

$\displaystyle\frac{1}{16}\operatorname*{\mathbb{E}}\left\|\nabla
g\right\|_{2}^{2}$
$\displaystyle\leq\sum_{r,s}((q-1)(q-2)B_{rs}^{(1)}+(q-1)B_{rs}^{(2)}+(q-1)B_{rs}^{(4)}+B_{rs}^{(5)})$
$\displaystyle\lesssim
q^{3}\operatorname{tr}^{2}(M)\left\|M\right\|_{F}^{2}+q^{2}\operatorname{tr}(M)\left\|M\right\|_{F}\sqrt{\operatorname{tr}(M^{4})}+q^{2}\left\|M\right\|_{F}^{4}+q\operatorname{tr}(M^{4}).$

By Poincaré’s inequality,

$\displaystyle\quad\ \operatorname*{Var}_{G}(\operatorname{tr}^{2}(G^{T}MG))$
(9) $\displaystyle\lesssim\operatorname*{\mathbb{E}}\left\|\nabla
g\right\|_{2}^{2}$ $\displaystyle\lesssim
q^{3}\operatorname{tr}^{2}(M)\left\|M\right\|_{F}^{2}+q^{2}\operatorname{tr}(M)\left\|M\right\|_{F}\sqrt{\operatorname{tr}(M^{4})}+q^{2}\left\|M\right\|_{F}^{4}+q\operatorname{tr}(M^{4}).$

Similar to before, each term on the right-hand side grows geometrically.

#### Step 2b.

Next we deal with
$\operatorname*{Var}_{M}(\operatorname*{\mathbb{E}}_{G}\operatorname{tr}^{2}\left(G^{T}MG\right)|M)$.

$\operatorname*{\mathbb{E}}\operatorname{tr}^{2}\left(G^{T}MG\right)=\operatorname*{\mathbb{E}}\bigg{(}\sum_{i,k,l}M_{kl}G_{ki}G_{li}\bigg{)}^{2}=\sum_{i,j}\sum_{k,l,k^{\prime},l^{\prime}}M_{kl}M_{k^{\prime}l^{\prime}}\operatorname*{\mathbb{E}}G_{ki}G_{li}G_{k^{\prime}j}G_{l^{\prime}j}.$

When $i\neq j$, for non-zero contribution, it must hold that $k=l$ and
$k^{\prime}=l^{\prime}$ and thus the nonzero contribution is

$\sum_{i\neq
j}\sum_{k,k^{\prime}}M_{kk}M_{k^{\prime}k^{\prime}}=q(q-1)\operatorname{tr}^{2}(M).$

When $i=j$, the contribution is (this is exactly the same as (7) in Step 1b.)

$\sum_{i}\sum_{\begin{subarray}{c}k,k^{\prime},l,l^{\prime}\end{subarray}}M_{kl}M_{k^{\prime}l^{\prime}}\operatorname*{\mathbb{E}}G_{ki}G_{li}G_{k^{\prime}i}G_{l^{\prime}i}=2q\left\|M\right\|_{F}^{2}+q\operatorname{tr}^{2}(M).$

Hence

$\operatorname*{\mathbb{E}}\operatorname{tr}^{2}\left(G^{T}MG\right)=2q\left\|M\right\|_{F}^{2}+q^{2}\operatorname{tr}^{2}(M)$

and when $M$ is random,

$\displaystyle\quad\
\operatorname*{Var}\left(\left.\operatorname*{\mathbb{E}}\operatorname{tr}^{2}(G^{T}MG)\right|M\right)$
(10)
$\displaystyle=\operatorname*{Var}\left(2q\left\|M\right\|_{F}^{2}+q^{2}\operatorname{tr}^{2}(M)\right)$
$\displaystyle\leq
4q^{2}\operatorname*{Var}(\left\|M\right\|_{F}^{2})+q^{4}\operatorname*{Var}(\operatorname{tr}^{2}(M))+2q^{3}\sqrt{\operatorname*{Var}(\left\|M\right\|_{F}^{2})\operatorname*{Var}(\operatorname{tr}^{2}(M))}.$

#### Step 3.

Let $U_{r}$ denote the variance of $\operatorname{tr}((A_{r}^{T}A_{r})^{2})$
and $V_{r}$ the variance of $\operatorname{tr}^{2}(A_{r}^{T}A_{r})$. Combining
(5), (6), (8), (9), (10), we have the following recurrence relations, where
$C_{1},C_{2},C_{3},C_{4}>0$ are absolute constants.

$\displaystyle U_{r+1}$ $\displaystyle\leq
C_{1}P_{r}+2U_{r}+\frac{1}{d_{r}^{2}}V_{r}+\frac{3}{d_{r}}\sqrt{U_{r}V_{r}}$
$\displaystyle V_{r+1}$ $\displaystyle\leq
C_{2}Q_{r}+\frac{1}{d_{r}^{2}}U_{r}+V_{r}+\frac{2}{d_{r}}\sqrt{U_{r}V_{r}}$
$\displaystyle P_{r+1}$ $\displaystyle\leq C_{3}P_{r}$ $\displaystyle Q_{r+1}$
$\displaystyle\leq C_{4}Q_{r}$ $\displaystyle U_{0}$ $\displaystyle=V_{0}=0$

In the base case, set $M=I_{p}$ (the $p\times p$ identity matrix in (6)) and
note that the second term in (5) vanishes. We see that
$P_{1}\lesssim(p^{3}q+pq^{3})/d_{1}^{4}$ after proper normalization.
(Alternatively we can calculate this precisely, see Appendix C.) Similarly we
have $Q_{1}\lesssim p^{3}q^{3}/d_{1}^{4}$. Note that
$Q_{1}/d_{1}^{2}\lesssim(p^{3}q+pq^{3})/d_{1}^{4}$. Now, we can solve that

$U_{r+1}\leq C^{r}\frac{p^{3}q+pq^{3}}{d_{1}^{4}}$

for some absolute constant $C>0$.

## References

  * [1] Gernot Akemann and Jesper R Ipsen.  Recent exact and asymptotic results for products of independent random matrices.  Acta Physica Polonica B, pages 1747–1784, 2015.
  * [2] Ahmed El Alaoui and Michael W. Mahoney.  Fast randomized kernel ridge regression with statistical guarantees.  In Proceedings of the 28th International Conference on Neural Information Processing Systems - Volume 1, NIPS’15, page 775–783, Cambridge, MA, USA, 2015. MIT Press.
  * [3] Richard Bellman.  Limit theorems for non-commutative operations. I.  Duke Mathematical Journal, 21(3):491–500, 1954.
  * [4] Giancarlo Benettin.  Power-law behavior of lyapunov exponents in some conservative dynamical systems.  Physica D: Nonlinear Phenomena, 13(1-2):211–220, 1984.
  * [5] Stéphane Boucheron, Gábor Lugosi, and Pascal Massart.  Concentration Inequalities: A Nonasymptotic Theory of Independence.  Oxford University Press, 2013.
  * [6] Matthew Brennan, Guy Bresler, and Dheeraj Nagaraj.  Phase transitions for detecting latent geometry in random graphs.  Probability Theory and Related Fields, pages 1215–1289, 2020.  doi:10.1007/s00440-020-00998-3.
  * [7] Sébastien Bubeck, Jian Ding, Ronen Eldan, and Miklós Z. Rácz.  Testing for high-dimensional geometry in random graphs.  Random Structures & Algorithms, 49(3):503–532, 2016.  doi:10.1002/rsa.20633.
  * [8] Sébastien Bubeck and Shirshendu Ganguly.  Entropic CLT and Phase Transition in High-dimensional Wishart Matrices.  International Mathematics Research Notices, 2018(2):588–606, 12 2016.  doi:10.1093/imrn/rnw243.
  * [9] Z. Burda, R. A. Janik, and B. Waclaw.  Spectrum of the product of independent random gaussian matrices.  Physical Review E, 81(4), Apr 2010.  doi:10.1103/physreve.81.041132.
  * [10] Emmanuel J Candes.  The restricted isometry property and its implications for compressed sensing.  Comptes rendus mathematique, 346(9-10):589–592, 2008.
  * [11] Andrea Crisanti, Giovanni Paladin, and Angelo Vulpiani.  Products of random matrices: in Statistical Physics, volume 104\.  Springer Science & Business Media, 2012.
  * [12] David L. Donoho.  Compressed sensing.  IEEE Transactions on Information Theory, 52(4):1289–1306, 2006\.
  * [13] S. Iida, H.A. Weidenmüller, and J.A. Zuk.  Statistical scattering theory, the supersymmetry method and universal conductance fluctuations.  Annals of Physics, 200(2):219–270, 1990.
  * [14] Jesper R. Ipsen.  Products of independent Gaussian random matrices.  PhD thesis, Bielefeld University, 2015.
  * [15] Hiroshi Ishitani.  A central limit theorem for the subadditive process and its application to products of random matrices.  Publications of the Research Institute for Mathematical Sciences, 12(3):565–575, 1977.
  * [16] Tiefeng Jiang.  How many entries of a typical orthogonal matrix can be approximated by independent normals?  The Annals of Probability, 34(4):1497–1529, Jul 2006.  doi:10.1214/009117906000000205.
  * [17] Tiefeng Jiang and Yutao Ma.  Distances between random orthogonal matrices and independent normals.  Transactions of the American Mathematical Society, 372(3):1509–1553, August 2019.  doi:10.1090/tran/7470.
  * [18] Ravindran Kannan and Santosh Vempala.  Spectral algorithms.  Now Publishers Inc, 2009.
  * [19] Michael Kapralov, Vamsi Potluru, and David Woodruff.  How to fake multiply by a gaussian matrix.  In International Conference on Machine Learning, pages 2101–2110. PMLR, 2016.
  * [20] Michael W. Mahoney.  Randomized algorithms for matrices and data.  Found. Trends Mach. Learn., 3(2):123–224, February 2011.  doi:10.1561/2200000035.
  * [21] Satya N Majumdar and Grégory Schehr.  Top eigenvalue of a random matrix: large deviations and third order phase transition.  Journal of Statistical Mechanics: Theory and Experiment, 2014(1):P01012, 2014.
  * [22] Robert M May.  Will a large complex system be stable?  Nature, 238(5364):413–414, 1972.
  * [23] P.A. Mello, P. Pereyra, and N. Kumar.  Macroscopic approach to multichannel disordered conductors.  Annals of Physics, 181(2):290–317, 1988.
  * [24] Ralf R Muller.  On the asymptotic eigenvalue distribution of concatenated vector-valued fading channels.  IEEE Transactions on Information Theory, 48(7):2086–2091, 2002\.
  * [25] Shanmugavelayutham Muthukrishnan.  Data streams: Algorithms and applications.  Now Publishers Inc, 2005.
  * [26] Guillaume Obozinski, Martin J Wainwright, Michael I Jordan, et al.  Support union recovery in high-dimensional multivariate regression.  The Annals of Statistics, 39(1):1–47, 2011.
  * [27] James C Osborn.  Universal results from an alternate random-matrix model for QCD with a baryon chemical potential.  Physical review letters, 93(22):222001, 2004.
  * [28] G Paladin and A Vulpiani.  Scaling law and asymptotic distribution of lyapunov exponents in conservative dynamical systems with many degrees of freedom.  Journal of Physics A: Mathematical and General, 19(10):1881, 1986\.
  * [29] Saurabh Paul, Christos Boutsidis, Malik Magdon-Ismail, and Petros Drineas.  Random projections for linear support vector machines.  ACM Transactions on Knowledge Discovery from Data (TKDD), 8(4):1–25, 2014.
  * [30] Miklós Z. Rácz and Jacob Richey.  A smooth transition from wishart to goe.  Journal of Theoretical Probability, pages 898–906, 2019.  doi:10.1007/s10959-018-0808-2.
  * [31] Garvesh Raskutti and Michael W Mahoney.  A statistical perspective on randomized sketching for ordinary least-squares.  The Journal of Machine Learning Research, 17(1):7508–7538, 2016\.
  * [32] Tamas Sarlos.  Improved approximation algorithms for large matrices via random projections.  In 2006 47th Annual IEEE Symposium on Foundations of Computer Science (FOCS’06), pages 143–152. IEEE, 2006.
  * [33] Khalid Shebrawi and Hussein Albadawi.  Trace inequalities for matrices.  Bulletin of the Australian Mathematical Society, 87(1):139–148, 2013.  doi:10.1017/S0004972712000627.
  * [34] Terence Tao.  Topics in random matrix theory, volume 132.  American Mathematical Soc., 2012.
  * [35] Roman Vershynin.  Introduction to the non-asymptotic analysis of random matrices.  In Yonina C. Eldar and Gitta Kutyniok, editors, Compressed Sensing: Theory and Applications, pages 210–268. Cambridge University Press, 2012.  doi:10.1017/CBO9780511794308.006.
  * [36] David P. Woodruff.  Sketching as a tool for numerical linear algebra.  Found. Trends Theor. Comput. Sci., 10(1–2):1–157, October 2014\.  doi:10.1561/0400000060.
  * [37] Fan Yang, Sifan Liu, Edgar Dobriban, and David P Woodruff.  How to reduce dimension with PCA and random projections?  arXiv:2005.00511 [math.ST], 2020.

##  Appendix A Omitted Calculations in Section 4.1

$\displaystyle
S_{1}(p,q)=3\sum_{i}\operatorname*{\mathbb{E}}B_{1i}^{4}+3\sum_{i\neq
j}\operatorname*{\mathbb{E}}B_{1i}^{2}B_{1j}^{2}=3dT_{1}(p,d)+3d(d-1)T_{4}(p,d)$
$\displaystyle\begin{aligned}
S_{3}(p,q)=\operatorname*{\mathbb{E}}A_{11}^{2}A_{21}^{2}&=\operatorname*{\mathbb{E}}\bigg{(}\sum_{i}B_{1i}G_{i1}\bigg{)}^{2}\bigg{(}\sum_{k}B_{2k}G_{k1}\bigg{)}^{2}\\\
&=\sum_{i,j,k,l}\operatorname*{\mathbb{E}}B_{1i}B_{1j}B_{2k}B_{2l}\operatorname*{\mathbb{E}}G_{i1}G_{j1}G_{k1}G_{l1}\\\
&=3\sum_{i}\operatorname*{\mathbb{E}}B_{1i}^{2}B_{2i}^{2}+\sum_{i\neq
j}\operatorname*{\mathbb{E}}B_{1i}^{2}B_{2j}^{2}+2\sum_{i\neq
j}\operatorname*{\mathbb{E}}B_{1i}B_{2i}B_{1j}B_{2j}\\\
&=3dT_{3}(p,d)+d(d-1)T_{5}(p,d)+2d(d-1)T_{6}(p,d)\end{aligned}$
$\displaystyle\begin{aligned}
S_{4}(p,q)=\operatorname*{\mathbb{E}}A_{11}^{2}A_{12}^{2}&=\operatorname*{\mathbb{E}}\bigg{(}\sum_{i}B_{1i}G_{i1}\bigg{)}^{2}\bigg{(}\sum_{k}B_{1k}G_{k2}\bigg{)}^{2}\\\
&=\sum_{i,j,k,l}\operatorname*{\mathbb{E}}B_{1i}B_{1j}B_{1k}B_{1l}\operatorname*{\mathbb{E}}G_{i1}G_{j1}G_{k2}G_{l2}\\\
&=\sum_{i}\operatorname*{\mathbb{E}}B_{1i}^{4}+\sum_{i\neq
j}\operatorname*{\mathbb{E}}B_{1i}^{2}B_{1j}^{2}=dT_{1}(p,d)+d(d-1)T_{4}(p,d).\end{aligned}$
$\displaystyle\begin{aligned}
S_{5}(p,q)=\operatorname*{\mathbb{E}}A_{11}^{2}A_{22}^{2}&=\operatorname*{\mathbb{E}}\bigg{(}\sum_{i}B_{1i}G_{i1}\bigg{)}^{2}\bigg{(}\sum_{k}B_{2k}G_{k2}\bigg{)}^{2}\\\
&=\sum_{i,j,k,l}\operatorname*{\mathbb{E}}B_{1i}B_{1j}B_{2k}B_{2l}\operatorname*{\mathbb{E}}G_{i1}G_{j1}G_{k2}G_{l2}\\\
&=\sum_{i,j}B_{1i}^{2}B_{2j}^{2}=dT_{3}(p,d)+d(d-1)T_{5}(p,d)\end{aligned}$
$\displaystyle\begin{aligned}
S_{6}(p,q)=\operatorname*{\mathbb{E}}A_{11}A_{12}A_{21}A_{22}&=\sum_{i,j,k,l}\operatorname*{\mathbb{E}}B_{1i}B_{1j}B_{2k}B_{2l}\operatorname*{\mathbb{E}}G_{i1}G_{j2}G_{k1}G_{l2}\\\
&=\sum_{i}\operatorname*{\mathbb{E}}B_{1i}^{2}B_{2i}^{2}+\sum_{i\neq
j}\operatorname*{\mathbb{E}}B_{1i}B_{1j}B_{2i}B_{2j}\\\
&=dT_{3}(p,d)+d(d-1)T_{6}(p,d)\end{aligned}$

##  Appendix B Omitted Calculations in Section 4.2

#### In Step 1a.

$\displaystyle\begin{aligned}
B^{(2)}_{r,s}&=\sum_{l,l^{\prime}}\sum_{u,v,v^{\prime}}M_{uv}M_{uv^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{us}^{2}X_{vj}X_{lj}X_{v^{\prime}j}X_{l^{\prime}j}\\\
&=\underbrace{\sum_{l\neq
l^{\prime}}\sum_{u}M_{ul}M_{ul^{\prime}}M_{rl}M_{rl^{\prime}}}_{v=l\neq
v^{\prime}=l^{\prime}}+\underbrace{\sum_{l}\sum_{\begin{subarray}{c}u\\\ v\neq
l\end{subarray}}M_{uv}^{2}M_{rl}^{2}}_{v=v^{\prime}\neq
l=l^{\prime}}+\underbrace{\sum_{l\neq
l^{\prime}}\sum_{u}M_{ul^{\prime}}M_{ul}M_{rl}M_{rl^{\prime}}}_{v=l^{\prime}\neq
l=v^{\prime}}\\\
&\qquad\qquad\qquad+3\underbrace{\sum_{l,u}M_{ul}^{2}M_{rl}^{2}}_{v=v^{\prime}=l=l^{\prime}}\\\
&=\left(\sum_{u,v}M_{uv}^{2}\right)\left(\sum_{l}M_{rl}^{2}\right)+2\sum_{l,l^{\prime},u}M_{ul^{\prime}}M_{ul}M_{rl}M_{rl^{\prime}}\\\
&=\left\|M\right\|_{F}^{2}\left\|M_{r,\cdot}\right\|_{2}^{2}+2\sum_{u}\langle
M_{u,\cdot},M_{r,\cdot}\rangle^{2}.\end{aligned}$
$\displaystyle\begin{aligned} B^{(3)}_{r,s}&=\sum_{j^{\prime}\neq
s}\left[\sum_{l,l^{\prime}}\sum_{u,v}M_{uv}M_{u^{\prime}l^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{us}X_{u^{\prime}s}X_{vs}X_{ls}X_{l^{\prime}j^{\prime}}^{2}\right]\\\
&=\underbrace{\sum_{l,l^{\prime}}\sum_{u\neq
l}M_{ul}M_{ul^{\prime}}M_{rl}M_{rl^{\prime}}}_{u=u^{\prime}\neq
v=l}+\underbrace{\sum_{l,l^{\prime}}\sum_{u\neq
l}M_{uu}M_{ll^{\prime}}M_{rl}M_{rl^{\prime}}}_{u=v\neq
u^{\prime}=l}+\underbrace{\sum_{l,l^{\prime}}\sum_{v}M_{lv}M_{vl^{\prime}}M_{rl}M_{rl^{\prime}}}_{u=l\neq
u^{\prime}=v}\\\
&\qquad+3\underbrace{\sum_{l,l^{\prime}}M_{ll}M_{ll^{\prime}}M_{rl}M_{rl^{\prime}}}_{u=u^{\prime}=v=l}\\\
&=\sum_{l,l^{\prime}}\sum_{u}(2M_{ul}M_{ul^{\prime}}+M_{uu}M_{ll^{\prime}})M_{rl}M_{rl^{\prime}}\\\
&=\sum_{l,l^{\prime}}(2\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle+\operatorname{tr}(M)M_{ll^{\prime}})M_{rl}M_{rl^{\prime}}.\end{aligned}$
$\displaystyle\begin{aligned}
\sum_{r,s}B^{(3)}_{r,s}&=\sum_{r,s}\sum_{l,l^{\prime}}(2\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle+\operatorname{tr}(M)M_{l^{\prime}l^{\prime}})M_{rl}M_{rl^{\prime}}\\\
&=q\sum_{l,l^{\prime}}(2\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle+\operatorname{tr}(M)M_{ll^{\prime}})\sum_{r}M_{rl}M_{rl^{\prime}}\\\
&=q\sum_{l,l^{\prime}}(2\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle+\operatorname{tr}(M)M_{ll^{\prime}})\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle\\\ &=q\sum_{l,l^{\prime}}2\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle^{2}+q\operatorname{tr}(M)\sum_{l,l^{\prime}}M_{ll^{\prime}}\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle\\\ &\leq
2q\operatorname{tr}(M^{4})+q\operatorname{tr}(M)\bigg{(}\sum_{l,l^{\prime}}M_{ll^{\prime}}^{2}\bigg{)}^{\frac{1}{2}}\bigg{(}\sum_{l,l^{\prime}}\langle
M_{l,\cdot},M_{l^{\prime},\cdot}\rangle^{2}\bigg{)}^{\frac{1}{2}}\\\ &\leq
2q\operatorname{tr}(M^{4})+q\operatorname{tr}(M)\left\|M\right\|_{F}\sqrt{\operatorname{tr}(M^{4})}\end{aligned}$

#### In Step 1b.

$\displaystyle\quad\ \sum_{i}\sum_{\begin{subarray}{c}k,l\\\
k^{\prime},l^{\prime}\end{subarray}}M_{kl}M_{k^{\prime}l^{\prime}}\operatorname*{\mathbb{E}}G_{ki}G_{li}G_{k^{\prime}i}G_{l^{\prime}i}$
$\displaystyle=\sum_{i}\bigg{(}\underbrace{\sum_{k\neq
l}M_{kl}^{2}}_{k=k^{\prime}\neq l=l^{\prime}}+\underbrace{\sum_{k\neq
l}M_{kl}^{2}}_{k=l^{\prime}\neq k^{\prime}=l}+\underbrace{\sum_{k\neq
l}M_{kk}M_{k^{\prime}k^{\prime}}}_{k=l\neq
k^{\prime}=l^{\prime}}+3\underbrace{\sum_{k}M_{kk}^{2}}_{k=k^{\prime}=l=l^{\prime}}\bigg{)}$
$\displaystyle=\sum_{i}\bigg{(}\sum_{k,l}M_{kl}^{2}+\sum_{k,l}M_{kl}^{2}+\sum_{k,l}M_{kk}M_{k^{\prime}k^{\prime}}\bigg{)}$
$\displaystyle=\sum_{i}(2\left\|M\right\|_{F}^{2}+\operatorname{tr}^{2}(M))$
$\displaystyle=2q\left\|M\right\|_{F}^{2}+q\operatorname{tr}^{2}(M).$

#### In Step 2a.

$\displaystyle B^{(2)}_{rs}$
$\displaystyle=\sum_{l,l^{\prime},u,u^{\prime},v}M_{uv}M_{u^{\prime}u^{\prime}}M_{rl}M_{rl^{\prime}}\operatorname*{\mathbb{E}}X_{ls}X_{l^{\prime}s}X_{us}X_{vs}X_{u^{\prime}j^{\prime}}^{2}$
$\displaystyle=\underbrace{\sum_{\begin{subarray}{c}l\neq u\\\
u^{\prime}\end{subarray}}M_{uu}M_{u^{\prime}u^{\prime}}M_{rl}^{2}}_{l=l^{\prime}\neq
u=v}+\underbrace{\sum_{\begin{subarray}{c}l\neq l^{\prime}\\\
u^{\prime}\end{subarray}}M_{ll^{\prime}}M_{u^{\prime}u^{\prime}}M_{rl}M_{rl^{\prime}}}_{l=u\neq
l^{\prime}=v}+\underbrace{\sum_{\begin{subarray}{c}l\neq l^{\prime}\\\
u^{\prime}\end{subarray}}M_{ll^{\prime}}M_{u^{\prime}u^{\prime}}M_{rl}M_{rl^{\prime}}}_{l=v\neq
l^{\prime}=u}$
$\displaystyle\qquad\qquad\qquad+3\underbrace{\sum_{l,u^{\prime}}M_{ll}M_{u^{\prime}u^{\prime}}M_{rl}^{2}}_{l=u=l^{\prime}=v}$
$\displaystyle=\operatorname{tr}(M)\left(\sum_{l,u}M_{uu}M_{rl}^{2}+2\sum_{l,l^{\prime}}M_{ll^{\prime}}M_{rl}M_{rl^{\prime}}\right)$
$\displaystyle=\operatorname{tr}(M)\left(\operatorname{tr}(M)\left\|M_{r,\cdot}\right\|_{2}^{2}+2\sum_{l,l^{\prime}}M_{ll^{\prime}}M_{rl}M_{rl^{\prime}}\right)$

$\displaystyle B^{(4)}_{rs}$
$\displaystyle=\sum_{l,u,u^{\prime},v,v^{\prime}}M_{uv}M_{u^{\prime}v^{\prime}}M_{rl}^{2}\operatorname*{\mathbb{E}}X_{uj}X_{vj}X_{u^{\prime}j}X_{v^{\prime}j}$
$\displaystyle=\left(\sum_{l}M_{rl}^{2}\right)\left(\underbrace{\sum_{u,v}M_{uv}^{2}}_{u=u^{\prime}\neq
v=v^{\prime}}+\underbrace{\sum_{u,v}M_{uv}^{2}}_{u=v^{\prime}\neq
u=v}+\underbrace{\sum_{u,v}M_{uv}^{2}}_{u=v\neq
u^{\prime}=v^{\prime}}+\underbrace{3\sum_{u}M_{uu}^{2}}_{u=v=u^{\prime}=v^{\prime}}\right)$
$\displaystyle=3\left(\sum_{l}M_{rl}^{2}\right)\sum_{u,v}M_{u,v}^{2}$
$\displaystyle=3\left\|M_{r,\cdot}\right\|_{2}^{2}\left\|M\right\|_{F}^{2}$

##  Appendix C Exact Variance when $r=2$

Suppose that $A$ is rotationally invariant under both left- and right-
multiplication of an orthogonal matrix. Define

$\displaystyle U_{1}(p,q)$
$\displaystyle=\operatorname*{Var}((A^{T}A)_{ii}^{2})$ $\displaystyle
U_{2}(p,q)$ $\displaystyle=\operatorname*{Var}((A^{T}A)_{ij}^{2})\quad i\neq
j$ $\displaystyle U_{3}(p,q)$
$\displaystyle=\operatorname{cov}((A^{T}A)_{ii}^{2},(A^{T}A)_{ik}^{2})\quad
i\neq k\quad\text{(same row, one entry on diagonal)}$ $\displaystyle
U_{4}(p,q)$
$\displaystyle=\operatorname{cov}((A^{T}A)_{ij}^{2},(A^{T}A)_{ik}^{2})\quad
j\neq k\quad\text{(same row, both entries off-diagonal)}$ $\displaystyle
U_{5}(p,q)$
$\displaystyle=\operatorname{cov}((A^{T}A)_{ii}^{2},(A^{T}A)_{jj}^{2})\quad
i\neq j\quad\text{(diff.\ rows and cols, both entries on diagonal)}$
$\displaystyle U_{6}(p,q)$
$\displaystyle=\operatorname{cov}((A^{T}A)_{ii}^{2},(A^{T}A)_{jk}^{2})\quad
i\neq j\neq k\quad\text{(diff.\ rows and cols, one entry on diagonal)}$
$\displaystyle U_{7}(p,q)$
$\displaystyle=\operatorname{cov}((A^{T}A)_{ij}^{2},(A^{T}A)_{kl}^{2})\quad
i\neq j\neq k\neq l\quad\text{(diff.\ rows and cols, nonsymmetric around
diag.)}$

It is clear that they are well-defined.

$\displaystyle\quad\ \operatorname*{Var}(\operatorname{tr}((A^{T}A)^{2}))$
$\displaystyle=\operatorname*{Var}\bigg{(}\sum_{i,j}(A^{T}A)_{ij}^{2}\bigg{)}$
$\displaystyle=\sum_{i,j,k,l}\operatorname{cov}((A^{T}A)_{ij}^{2},(A^{T}A)_{kl}^{2})$
$\displaystyle=\sum_{i,j}\operatorname*{Var}((A^{T}A)_{ij}^{2})+2\sum_{i}\sum_{j\neq
l}\operatorname{cov}((A^{T}A)_{ij}^{2},(A^{T}A)_{il}^{2})+\sum_{\begin{subarray}{c}i\neq
k\\\ j\neq
l\end{subarray}}\operatorname{cov}(\operatorname*{\mathbb{E}}(A^{T}A)_{ij}^{2},(A^{T}A)_{kl}^{2})$
$\displaystyle=q\operatorname*{Var}((A^{T}A)_{11}^{2})+q(q-1)\operatorname*{Var}(\operatorname*{\mathbb{E}}(A^{T}A)_{12}^{2})$
$\displaystyle\qquad+2\left[2q(q-1)\operatorname{cov}((A^{T}A)_{11}^{2},(A^{T}A)_{12}^{2})+q(q-1)(q-2)\operatorname{cov}((A^{T}A)_{12}^{2},(A^{T}A)_{13}^{2})\right]$
$\displaystyle\qquad+q(q-1)\operatorname{cov}(\operatorname*{\mathbb{E}}(A^{T}A)_{11}^{2},(A^{T}A)_{22}^{2})+q(q-1)\operatorname{cov}(\operatorname*{\mathbb{E}}(A^{T}A)_{12}^{2},(A^{T}A)_{21}^{2})$
$\displaystyle\qquad\quad+2q(q-1)(q-2)\operatorname{cov}((A^{T}A)_{11}^{2},(A^{T}A)_{23}^{2})$
$\displaystyle\qquad\quad+2q(q-1)(q-2)\operatorname{cov}((A^{T}A)_{12}^{2},(A^{T}A)_{31}^{2})$
$\displaystyle\qquad\quad+q(q-1)(q-2)(q-3)\operatorname*{\mathbb{E}}(A^{T}A)_{12}^{2}(A^{T}A)_{34}^{2}$
$\displaystyle=qU_{1}(p,q)+q(q-1)U_{2}(p,q)+2q(q-1)(2U_{3}(p,q)+(q-2)U_{4}(p,q))$
$\displaystyle\qquad+q(q-1)(U_{5}(p,q)+U_{2}(p,q))+2q(q-1)(q-2)(U_{6}(p,q)+U_{4}(p,q))$
$\displaystyle\qquad+q(q-1)(q-2)(q-3)U_{7}(p,q)$
$\displaystyle=qU_{1}(p,q)+q(q-1)(2U_{2}(p,q)+4U_{3}(p,q)+U_{5}(p,q))$
$\displaystyle\qquad+2q(q-1)(q-2)(2U_{4}(p,q)+U_{6}(p,q))+q(q-1)(q-2)(q-3)U_{7}(p,q).$

Let us calculate $U_{1},\dots,U_{7}$ for a $p\times q$ Gaussian random matrix
$G$.

$\displaystyle\begin{aligned}
U_{1}(p,q)=\operatorname*{\mathbb{E}}(G^{T}G)_{11}^{4}-(\operatorname*{\mathbb{E}}(G^{T}G)_{11}^{2})^{2}&=\operatorname*{\mathbb{E}}\left\|G_{1}\right\|_{2}^{8}-(\operatorname*{\mathbb{E}}\left\|G_{1}\right\|_{2}^{4})^{2}\\\
&=p(p+2)(p+4)(p+6)-(p(p+2))^{2}\\\ &=8p(p+2)(p+3)\end{aligned}$
$\displaystyle\begin{aligned}
U_{2}(p,q)=\operatorname*{\mathbb{E}}(G^{T}G)_{12}^{4}-(\operatorname*{\mathbb{E}}(G^{T}G)_{12}^{2})^{2}&=\operatorname*{\mathbb{E}}\left(\sum_{r}G_{r1}G_{r2}\right)^{4}-(\operatorname*{\mathbb{E}}\langle
G_{1},G_{2}\rangle^{2})^{2}\\\
&=\sum_{r,s,t,u}\operatorname*{\mathbb{E}}G_{r1}G_{s1}G_{t1}G_{u1}G_{r2}G_{s2}G_{t2}G_{u2}-p^{2}\\\
&=3\sum_{r\neq
t}\operatorname*{\mathbb{E}}G_{r1}^{2}G_{t1}^{2}G_{r2}^{2}G_{t2}^{2}+\sum_{r}G_{r1}^{4}G_{r2}^{4}-p^{2}\\\
&=3p(p-1)+9p-p^{2}=2p(p+3).\end{aligned}$ $\displaystyle\begin{aligned}
U_{3}(p,q)&=\operatorname*{\mathbb{E}}(G^{T}G)_{11}^{2}(G^{T}G)_{12}^{2}-\operatorname*{\mathbb{E}}(G^{T}G)_{11}^{2}\operatorname*{\mathbb{E}}(G^{T}G)_{12}^{2}\\\
&=\operatorname*{\mathbb{E}}(G_{1}^{T}G_{1})^{2}G_{1}^{T}G_{2}G_{2}^{T}G_{1}-\operatorname*{\mathbb{E}}\left\|G_{1}\right\|_{2}^{4}\operatorname*{\mathbb{E}}\langle
G_{1},G_{2}\rangle^{2}\\\
&=\operatorname*{\mathbb{E}}(G_{1}^{T}G_{1})^{2}G_{1}^{T}(\operatorname*{\mathbb{E}}G_{2}G_{2}^{T})G_{1}-p(p+2)\cdot
p\\\ &=\operatorname*{\mathbb{E}}(G_{1}^{T}G_{1})^{3}-p^{2}(p+2)\\\
&=\operatorname*{\mathbb{E}}\left\|G_{1}\right\|_{2}^{6}-p^{2}(p+2)=p(p+2)(p+4)-p^{2}(p+2)=4p(p+2)\end{aligned}$
$\displaystyle\begin{aligned}
U_{4}(p,q)&=\operatorname*{\mathbb{E}}(G^{T}G)_{12}^{2}(G^{T}G)_{13}^{2}-\operatorname*{\mathbb{E}}(G^{T}G)_{12}^{2}\operatorname*{\mathbb{E}}(G^{T}G)_{13}^{2}\\\
&=\operatorname*{\mathbb{E}}G_{1}^{T}G_{2}G_{2}^{T}G_{1}G_{1}^{T}G_{3}G_{3}^{T}G_{1}-p^{2}\\\
&=\operatorname*{\mathbb{E}}G_{1}^{T}\operatorname*{\mathbb{E}}(G_{2}G_{2}^{T})G_{1}G_{1}^{T}\operatorname*{\mathbb{E}}(G_{3}G_{3}^{T})G_{1}-p^{2}\\\
&=\operatorname*{\mathbb{E}}(G_{1}^{T}G_{1})^{2}-p^{2}=\operatorname*{\mathbb{E}}\left\|G_{1}\right\|_{2}^{4}-p^{2}=p(p+2)-p^{2}=2p\end{aligned}$
$\displaystyle\begin{aligned}
U_{5}(p,q)&=\operatorname*{\mathbb{E}}(G^{T}G)_{11}^{2}(G^{T}G)_{22}^{2}-\operatorname*{\mathbb{E}}(G^{T}G)_{11}^{2}\operatorname*{\mathbb{E}}(G^{T}G)_{22}^{2}\\\
&=\operatorname*{\mathbb{E}}\left\|G_{1}\right\|_{2}^{4}\left\|G_{2}\right\|_{2}^{4}-\operatorname*{\mathbb{E}}\left\|G_{1}\right\|_{2}^{4}\left\|G_{2}\right\|_{2}^{4}=0\end{aligned}$
$\displaystyle\begin{aligned}
U_{6}(p,q)&=\operatorname*{\mathbb{E}}(G^{T}G)_{11}^{2}(G^{T}G)_{23}^{2}-\operatorname*{\mathbb{E}}(G^{T}G)_{11}^{2}\operatorname*{\mathbb{E}}(G^{T}G)_{23}^{2}\\\
&=\operatorname*{\mathbb{E}}\left\|G_{1}\right\|_{2}^{4}\langle
G_{2},G_{3}\rangle^{2}-\operatorname*{\mathbb{E}}\left\|G_{1}\right\|_{2}^{4}\operatorname*{\mathbb{E}}\langle
G_{2},G_{3}\rangle^{2}=0\end{aligned}$ $\displaystyle\begin{aligned}
U_{7}(p,q)&=\operatorname*{\mathbb{E}}(G^{T}G)_{12}^{2}(G^{T}G)_{34}^{2}-\operatorname*{\mathbb{E}}(G^{T}G)_{12}^{2}\operatorname*{\mathbb{E}}(G^{T}G)_{34}^{2}\\\
&=\operatorname*{\mathbb{E}}\langle G_{1},G_{2}\rangle^{2}\langle
G_{3},G_{4}\rangle^{2}-\operatorname*{\mathbb{E}}\langle
G_{1},G_{2}\rangle^{2}\operatorname*{\mathbb{E}}\langle
G_{3},G_{4}\rangle^{2}=0\end{aligned}$

Therefore

$\displaystyle\operatorname*{Var}(\operatorname{tr}((G^{T}G)^{2}))$
$\displaystyle=qU_{1}+q(q-1)(2U_{2}+4U_{3}+U_{5})+2q(q-1)(q-2)(2U_{4}+U_{6})$
$\displaystyle\qquad+q(q-1)(q-2)(q-3)U_{7}$
$\displaystyle=qU_{1}+q(q-1)(2U_{2}+4U_{3})+4q(q-1)(q-2)U_{4}$
$\displaystyle=4pq(5+5p+5q+2p^{2}+5pq+2q^{2}).$

When $r=2$, recalling that
$\operatorname*{\mathbb{E}}(A_{2}-A_{1})=(1+o(1))p^{2}q^{2}/d^{3}$ (see (4)),
we have that

$\frac{\sqrt{\operatorname*{Var}(\operatorname{tr}((\frac{1}{\sqrt{d}}G^{T}\cdot\frac{1}{\sqrt{d}}G)^{2}))}}{p^{2}q^{2}/d^{3}}\leq\frac{6d}{\max\\{p,q\\}^{\frac{1}{2}}\min\\{p,q\\}^{\frac{3}{2}}}.$

If the right-hand side above is at most a small constant $c$, we can
distinguish $A_{2}$ from $A_{1}$ with probability at least a constant.