The Product of Gaussian Matrices is Close to Gaussian

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.

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,

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.,

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,

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}$,

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

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

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.