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})$$

rowrank. Apologies! $\endgroup$ – student Jan 3 at 16:07orthonormalfor $(AB)^+=B^+ A^+$ to hold true. $\endgroup$ – Ludwig Jan 3 at 16:53