As an accidental byproduct of some numerical simulations I have been doing as part of a research paper in machine learning, I made the observation that the singular values of the random matrix $\frac{1}{T}XY^T$ are strikingly similar to the Marchenko-Pastur distribution with $\gamma=\frac{N}{2T}$ (Note the title didn't include the $\frac{1}{T}$ term for brevity, but it is necessary). Here $X,Y \in \mathbb{R}^{N \times T}$ with i.i.d. gaussian entries, with mean 0 and standard deviation $\sigma$ I.e. the two matrices X and Y are samples from the same exact ensemble.
Above is a plot showing MP(1/8) vs simulations of eigenvalues values of $\frac{1}{T^2}Y^TYX^TX$ and $\frac{1}{2T}[X|Y]^T[X|Y]$ that I did with N=250 and T=1000, averaged over 10 runs. The orange curve is the theoretical distribution for the eigenvalues of $\frac{1}{2T}AA^T, A \in \mathbb{R}^{N \times 2T}$, which we know is approximated by the Marchenko-Pastur distribution with $\gamma=\frac{N}{2T}=\frac{1}{8}$. The green series shows a simulation of these eigenvalues. It appears to match the theoretical limit arbitrarily well. The blue series is the simulated density of singular values of of $\frac{1}{T}XY^T,$ or the eigenvalues of $\frac{1}{T^2}YX^TXY^T$. I am not showing the plot of the theoretical density for that matrix here (as described for example in "Large dimension forecasting models and random singular value spectra" (2004) by Bouchaud et. al.) as it isn't relevant to my inquiry, but they do note that this is the same as the eigenvalues of $\frac{1}{T^2}Y^TYX^TX$. The main point is that, for any N<T for N,T sufficiently large, and for any $\sigma$ sufficiently small, the former distribution appears to very closely match the latter, if just slightly shifted to the left. Though it should be noted that, even when I control the seeds, there isn't a clear 1-1 correspondence between the eigenvalues of the two distributions (their empirical distributions are not literally related by a shift as you can see by the fact that the histograms are not parallel on a granular level -- rather their limiting distributions are closely parallel). This appears to hold even if I vary N and T, and $\sigma$.
My question is, why? Are they asymptotically equal in some limit? If not, and there is some persistent "error" term, can we identify the offending term, and at least say something about it's order?
In order to try to explain this, I wrote out an expression for the the singular values of $\frac{1}{\sqrt{2T}}[X|Y]$ which I will call M, noting that it dimension $(N \times 2T)$ and is identically distributed as $\frac{1}{\sqrt{2T}}A $. The singular values of M are the eigenvalues of the matrix $$M^TM=\frac{1}{2T}\begin{bmatrix} [X^TX] [Y^TX] \newline [X^TY] [Y^TY] \end{bmatrix} \in \mathbb{R}^{2T \times 2T},$$
which has characteristic equation $$0=Det(M^TM-I\lambda)=Det(\begin{bmatrix} [\frac{X^TX}{2T}-I \lambda] [\frac{Y^TX}{2T}] \newline [\frac{X^TY}{2T}] [\frac{Y^TY}{2T} - I \lambda] \end{bmatrix})$$
(A brief aside here: I could have equivalently looked at the singular values of M^T, or the eigenvalues of $MM^T \in \mathbb{R}^{Nx \times N} = \frac{1}{2T}XX^T + \frac{1}{2T}YY^T$. Perhaps this route could yield a result based on previous results on the spectra of free convolutions of free random matrices having additive R-Transforms, but it would likely be a more complex proof, so I chose the other route.)
Since $\frac{X^TX}{2T}-I \lambda$ is positive definite for $N<T$ (assuming a fixed $\lambda$ that is an eigenvalue of $M^TM$ is not also an eigenvalue of $\frac{X^TX}{2T}$), it is invertible and therefore (by a well-known identity), the characteristic equation can be rewritten as
$$Det(\frac{X^TX}{T}- I \lambda)Det(\frac{1}{T}Y^TY - I \lambda-\frac{X^TY}{T}(\frac{X^TX}{T}-I \lambda)^{-1}\frac{Y^TX}{T})=0$$
The next step I believe should be to make an argument that the term $\frac{X^TY}{T}(\frac{X^TX}{T}-I \lambda)^{-1}\frac{Y^TX}{T}$ is "small" relative to $\frac{1}{T}Y^TY - I \lambda$, since $||(\frac{1}{T}X^TX-I \lambda)^{-1}||$ and $||(\frac{1}{T}Y^TY-I \lambda)||$ are both O(1) (by Marchenko-Pastur), but $||\frac{1}{T}X^TY||,||\frac{1}{T}Y^TX|| \to O(\frac{\sigma^2 N}{T})$, which we get by calculating the variance of each entry in $X^TY$ and then applying an argument by Van Handel. Thus, as long as $\sigma^2<<\frac{T}{N}$, by Cauchy-Schwartz, $||\frac{1}{T}X^TY(\frac{1}{T}X^TX-I \lambda)^{-1}\frac{1}{T}Y^TX|| \to 0$
That seemed like a good start, but unfortunately I don't know where to go exactly from here. The above seems to tell us that in the limit of small $\sigma$, if $\lambda$ is an eigenvalue of either $\frac{1}{T}X^TX$ or $\frac{1}{T}Y^TY$, then it is also an eigenvalue of $M^TM$. But this isn't particularly useful. First of all, we are looking for a relationship between the eigenvalues of $M^TM$ and $\frac{1}{T}Y^TYX^TX$, not $\frac{1}{T}X^TX$ or $\frac{1}{T}Y^TY$ individually. Eliminating that small term by allowing $\sigma \to 0$ in simulations does not get rid of the difference, so I conclude the difference does not come from the approximation error related to $\sigma, N$ and $T$ described above.
I could multply out the expression as
$$Det(\frac{1}{T}X^TX- I \lambda)Det(\frac{1}{T}Y^TY - I \lambda) = Det(\frac{1}{T^2}X^TXY^TY - \frac{1}{T}\lambda(X^TX+Y^TY)+\lambda^2 I),$$
Which is tempting because the term $\frac{1}{T^2}X^TXY^TY$ appears inside the determinant. However, Is the goal here is to try to form some kind of bijection between individual eigenvalues of the two matrices, or rather to say something about the similarity of their distributions more generally? Proving that $\lambda$ is an eigenvalue of $M^TM \implies \exists! \lambda^*$ "close" to $\lambda$ s.t. $\lambda^*$ is an eigenvalue of $\frac{1}{T^2}Y^TYX^TX$ would imply the distributions are similar, but perhaps I can prove the latter without the former?
Any help on completing this proof would be greatly appreciated.
