Let $A\in\mathbb{R}^{r\times m}$ be a matrix of full row rank, and let $\cdot^+$ denote the Moore-Penrose inverse.
Consider a sequence of matrices $\{B_n\}_{n>1}$, $B_n\in\mathbb{R}^{m\times n}$, whose entries are i.i.d. random Gaussian variables with zero mean and finite variance. Numerical simulations suggest that $$ \|(AB_n)^+ - B_n^+ A^+\| \to 0\ \ \text{ as }\ \ n\to\infty $$ where $\|\cdot\|$ denotes the 2-norm of a matrix.
Is it possible to provide a formal proof of this claim?
My question could be either very silly (if so, please close it) or a well-known fact (if so, I will be glad if you can provide pointers to the literature). In any case, the observed numerical behavior looks quite surprising to me, since it is rather well-known that for $A$ and $B$ both of full row rank it holds $(AB)^+ \ne B^+ A^+$ (except infor some very special cases).