Not much literature seem to address this connection. Is it fair to say that the least square linear regression (commonly seen in machine learning) is an approximation of the linear MMSE estimator (well-known in estimation theory)?
(1) The LS linear regression solves for
$$\pmb{\hat\beta}_{LS} = \underset{\pmb\beta}{\textrm{argmin}} \: \| \mathbf{y-X\pmb\beta\|^2},$$
where $\mathbf y$ is a vector of $n$ outputs of the training set, and $\mathbf X$ is an $n\times p$ matrix consists of $n$ input vectors of the training set. Solving it, we have
$$\pmb{\hat \beta}_{LS} = \mathbf{(X^TX)^{-1}X^Ty}.$$
(2) The linear MMSE estimator, on the other hand, seeks
$$\pmb{\hat \beta}_{LMMSE} = \underset{\pmb\beta}{\textrm{argmin}} \: E\left[ (Y-\pmb{X^T\beta})^2\right],$$
where $Y$ is the output random variable that we want to estimate, and $\pmb X$ is the $p$-dimensional input random vector. This leads to
$$\pmb{\hat \beta}_{LMMSE} = \mathbf R_{\pmb X}^{-1}\mathbf R_{\pmb XY}.$$
(3) From (1) & (2), it's clear that if the training set of (1) is ergodic in the mean, then
$$\frac{1}{n} \mathbf{X^TX} \to \mathbf R_{\pmb X} \:\mathrm{(m.s.)}, $$ $$\frac{1}{n} \mathbf{X^Ty} \to \mathbf R_{\pmb XY} \:\mathrm{(m.s.),}$$
as $n\to \infty$. Moreover, since matrix inversion is a continuous mapping, we also have $n \mathbf{(X^TX)^{-1}} \to \mathbf R_{\pmb X}^{-1}.$ Therefore,
$$\pmb{\hat \beta}_{LS} = \mathbf{(X^TX)^{-1}X^Ty} \quad \to \quad \mathbf R_{\pmb X}^{-1}\mathbf R_{\pmb XY} = \pmb{\hat \beta}_{LMMSE}.$$
Was I mistaken somewhere above? Are there pitfalls or caveats that I should heed when making this connection? Thanks a lot!