Spectrum asymptotics for a product of $k$ random matrices?

Question

How does the spectrum of a product of $k$ random matrices behave around 0?

In particular, I'm wondering if the CDF of squared singular values behaves as $x^{\frac{1}{k+1}}$ around 0. The result for $k=1$ comes from Marchenko-Pastur.

Below is a visualization of empirical CDF of squared singular values of a matrix formed by multiplying together $k$ square $1000\times 1000$ matrices with standard Normal entries.

Notebook

Unclear from simulation above is whether apparent "faster than power-law" decay is a finite sample effect or if it persists for $n\to\infty$. Any tips?

Motivation: a potential candidate for RMT model to match scaling laws observed in neural network training (background). Shape of loss curve over time for $t\to \infty$ comes from shape of spectrum for $x\to 0$

Answer 1

The distribution of the singular values has been studied in Singular value correlation functions for products of Wishart random matrices and Eigenvalues and Singular Values of Products of Rectangular Gaussian Random Matrices. Near the origin the density of the squared singular values $x$ diverges as $\rho(x)\propto x^{-k/(k+1)}$, so a CDF that scales $\propto x^{1/(k+1)}$, as anticipated in the OP.

For completeness, I mention that the distibution of the eigenvalues of the product of $k$ random matrices has been studied in Spectrum of the Product of Independent Random Gaussian Matrices. It is circularly symmetric in the complex plane, with density $\rho(z)\propto |z|^{-2+2/k}$ in a disc of finite radius determined by the variance of the matrix elements. For $k=1$ one thus recovers the uniform density of the Ginibre ensemble.