$(AB)^+\approx B^+A^+$ for $B$ "fat" enough? 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 for some very special cases). 
 A: As @FedericoPoloni pointed out, this must hinge on the fact that the rows of $B_n$ tend to be orthogonal as $n$ increases. In fact, $$\mathrm{E}[(B_n B_n^*)_{ij}] = n \sigma^2 \delta_{ij} \\ \mathrm{Cov}[(B_n B_n^*)_{ij}, (B_n B_n^*)_{i'j'}] = n \sigma^4 \, (\delta_{ii'} \delta_{jj'} + \delta_{ij'} \delta_{i'j})$$ so that we might as well write $$B_n B_n^* = n \sigma^2 I + \sqrt{n} \sigma^2 R$$ where $R$ is a random matrix with entries of zero mean and $O(1)$ variance. The near-orthogonality comes into play as $$\lim_{n \to \infty} \frac{1}{n} B_n B_n^* = \sigma^2 I$$ which should feature somehow in showing the suggestion.
To introduce this into the M-P inverse, the only thing that comes to mind right now is to take the limit definition $$A^+ = \lim_{\delta \searrow 0} A^* (A A^* + \delta I)^{-1}$$ to define a sequence $$(A B_n)^+_k = B_n^* A^* (A B_n B_n^* A^* + \tfrac{1}{k} I)^{-1}$$ that converges to $(A B_n)^+$. We also have $B_n^+ = B_n^* (B_n B_n^*)^{-1}$ a.s. for $n \ge m$ and $A^+ = A^* (A A^*)^{-1}$. Hence for a fixed $k$ we can evaluate $$\Vert (A B_n)^+_k - B_n^+ A^+ \Vert = \frac{1}{n} \left\Vert B_n^* A^* \left(A \, (\tfrac{1}{n} B_n B_n^*) \, A^* + \tfrac{1}{n k} I\right)^{-1} - B_n^* (\tfrac{1}{n} B_n B_n^*)^{-1} A^* (A A^*)^{-1} \right\Vert$$ which looks harmless enough to be pushed below any desired bound with a suitably large $n$.
If that's the case, pick $\epsilon > 0$ and use the triangle inequality $$\Vert (A B_n)^+ - B_n^+ A^+ \Vert \le \Vert (A B_n)^+ - (A B_n)^+_k \Vert + \Vert (A B_n)^+_k - B_n^+ A^+ \Vert$$ to look at each term on the rhs. individually: first find $k$ s.t. the first term is less than $\epsilon/2$, and then find $n$ s.t. the second term is less than $\epsilon/2$.

Moments of the entries of $BB^*$
For iid. normal entries $B_{ij}$ with zero mean and variance $\sigma^2$: $$\mathrm{E}[B_{ij} B_{kl}] = \sigma^2 \delta_{ik} \delta_{jl}$$
Hence the expectation: $$\mathrm{E}[(B B^*)_{ij}] = \sum_k \mathrm{E}[B_{ik} B_{jk}] = \sum_k \sigma^2 \delta_{ij} \delta_{kk} = n \sigma^2 \delta_{ij}$$
For jointly normal $X_1, X_2, X_3, X_4$ with zero mean: $$\mathrm{E}[X_1 X_2 X_3 X_4] = \mathrm{E}[X_1 X_2] \mathrm{E}[X_3 X_4] + \mathrm{E}[X_1 X_3] \mathrm{E}[X_2 X_4] + \mathrm{E}[X_1 X_4] \mathrm{E}[X_2 X_3]$$
Hence the covariance: $$
\mathrm{Cov}[(B B^*)_{ij}, (B B^*)_{i'j'}]
= \sum_{k,k'} \mathrm{Cov}[B_{ik} B_{jk}, B_{i'k'} B_{j'k'}] \\
= \sum_{k,k'} \left\{\mathrm{E}[B_{ik} B_{jk} B_{i'k'} B_{j'k'}] - \mathrm{E}[B_{ik} B_{jk}] \mathrm{E}[B_{i'k'} B_{j'k'}]\right\} \\
= \sum_{k,k'} \left\{\mathrm{E}[B_{ik} B_{i'k'}] \mathrm{E}[B_{jk} B_{j'k'}] + \mathrm{E}[B_{ik} B_{j'k'}] \mathrm{E}[B_{jk} B_{i'k'}]\right\} \\
= \sum_{k,k'} \left\{\sigma^4 \delta_{ii'} \delta_{kk'} \delta_{jj'} \delta_{kk'} + \sigma^4 \delta_{ij'} \delta_{kk'} \delta_{i'j} \delta_{kk'}\right\} \\
= n \sigma^4 (\delta_{ii'} \delta_{jj'} + \delta_{ij'} \delta_{i'j})$$
