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.

MP(1/8) vs simulations of the two distributions in question

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.


1 Answer 1


The distribution function of the singular values was calculated in Eigenvalues and Singular Values of Products of Rectangular Gaussian Random Matrices, see eq. (61). The moment generating function $M(z)=\sum_{n\geq 1 }m_n/z^n$ with $m_n$ the $n$-th moment of the singular values of $XY^\top$ (with $X,Y\in\mathbb{R}^{N\times T}$) is given by $$\frac{1}{M(z)}\bigl(1+M(z)\bigr)^2\bigl(1+(N/T)M(z)\bigr)=z/\sigma^2.$$ Marchenko-Pastur (MP) for the singular values of $Z\in\mathbb{R}^{N\times 2T}$ gives $$\frac{\sqrt{2T/N}}{M_{\text{MP}}(z)}\bigl(1+M_{\text{MP}}(z)\bigr)\bigl(1+(2T/N)M_{\text{MP}}(z)\bigr)=z/\sigma^2.$$

For example, if $N/T=1/4$ and $\sigma=1$ I find $$M(z)=\frac{1}{z} +\frac{9}{4 z^2}+\frac{105}{16 z^3}+\frac{1393}{64 z^4}+{\cal O}(z^{-5}),$$ while $$M_{\text{MP}}(z)=\frac{2 \sqrt{2}}{z}+\frac{72}{z^2}+\frac{1424 \sqrt{2}}{z^3}+\frac{60480}{z^4}+{\cal O}(z^{-5}).$$ A rescaling of $z$ with a factor $\sqrt{2T/N}=\sqrt 8$ aligns the first moment, but not higher moments.
  • $\begingroup$ Interesting. So can we use this to show that the two distributions are close in some way? I was hoping that comparing the moment generating function of Marchenko-Pastur for a $Nx2T$ random matrix, given as you say by $$\frac{1}{M(z)}\bigl(1+M(z)\bigr)\bigl(1+(N/2T)M(z)\bigr)=z/\sigma^2.,$$ would illuminate why these distributions are similar. Alas, when both are written out as polynomials in M, none of the terms match. So I'm still not sure why they are close, or how this closeness can be quantified. $\endgroup$
    – Paul
    Commented Aug 15, 2023 at 17:49
  • $\begingroup$ ah, I see that I need to take the first moment and then compare ... will get back once I complete that $\endgroup$
    – Paul
    Commented Aug 15, 2023 at 21:33
  • $\begingroup$ so in order to get the expectation, I first need to solve for M(0), and then write M'(0) in terms of M(0)? $\endgroup$
    – Paul
    Commented Aug 15, 2023 at 22:07
  • $\begingroup$ My conjecture from this perspective would be that the first moment of one distribution is slightly left shifted from the other, but the other moments are all the same. $\endgroup$
    – Paul
    Commented Aug 15, 2023 at 22:39
  • $\begingroup$ Thank you for your clarifying edits, Carlo. I'm surprised that the conclusion is the exact opposite of what I thought! They seemed to clearly have different means in my simulations, but the same overall shape. I guess looks can be deceiving! $\endgroup$
    – Paul
    Commented Aug 16, 2023 at 15:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.