Why doesn't correlation of residuals matter when testing for normality? When $Y = AX + \varepsilon$ (i.e., $Y$ comes from linear regression model),
$$\varepsilon \sim \mathcal{N}(0, \sigma^2 I) \hspace{1em} \Rightarrow
\hspace{1em} \hat{e} = (I - H) Y \sim \mathcal{N}(0, (I - H) \sigma^2_{})$$
and in that case residuals $\hat{e}_1, \ldots, \hat{e}_n$ are correlated and not independent. But when we do regression diagnostics and want to test the assumption
$\varepsilon \sim \mathcal{N}(0, \sigma^2 I)$, every textbook suggests to use
Q–Q plots and statistical tests on residuals $\hat{e}$ that were designed to
test whether $\hat{e} \sim \mathcal{N}(0, \sigma^2 I)$ for some $\sigma^2 \in
\mathbb{R}$. 
How come it doesn't matter for these tests that residuals are correlated, and
not independent? It is often suggested to use standardised residuals:
$$\hat{e}_i' = \frac{\hat{e}_i}{\sqrt{1 - h_{ii}}},$$
but that only makes them homoscedastic, not independent. 
To rephrase the question: Residuals from OLS regression are correlated. I understand that in practice, these correlations are so small (most of the time? always?), they can be ignored when testing whether residuals came from normal distribution. My question is, why?
 A: In your notation, $H$ is the projection an the column space of $X$, i.e. the subspace spanned of all regressors. Therefore $M:=I_{n}-H$ is the projection on everything orthogonal to the subspace spanned by all regressors. 
If $X\in\mathbb{R}^{n\times k}$, then $\hat{e}\in\mathbb{R}^{n}$ is singular normal distributed and the elements are correlated, as you state. 
The errors $\varepsilon$ are unobservable and are in general not orthogonal to the subspace spanned by $X$. 
For the sake of argument, assume that the error $\varepsilon\perp\operatorname{span}\left(X\right)$.
If this was true, we would have $y=X\beta+\varepsilon=\tilde{y}+\varepsilon$ with $\tilde{y}\perp\varepsilon$. Since $\tilde{y}=X\beta\in\operatorname{span}\left(X\right)$, we could decompose $y$ and get the true $\varepsilon$. 
Assume we have a basis $b_{1},\ldots,b_{n}$ of $\mathbb{R}^{n}$, where the first $b_{1},\ldots,b_{k}$ basis vector span the subspace $\operatorname{span}\left(X\right)$ and the remaining $b_{k+1},\ldots,b_{n}$ span $\operatorname{span}\left(X\right)^{\perp}$.
In general, the error $\varepsilon=\alpha_{1}b_{1}+\ldots+\alpha_{n}b_{n}$ will have non-zero components $\alpha_{i}$ for $i\in\left\{1,\ldots,k\right\}$. This non-zero components will get mixed up with $X\beta$ and therefore can not be recovered by projection on $\operatorname{span}\left(X\right)$. 
Since we can never hope to recover the true errors $\varepsilon$ and $\hat{e}$ are correlated singular $n$-dimensional normal, we could transform $\hat{e}\in\mathbb{R}^{n}\mapsto e^{*}\in\mathbb{R}^{n-k}$. There we can have that
\begin{equation}
e^{*}\sim\mathcal{N}_{n-k}\left(0,\sigma^{2}I_{n-k}\right)
\textrm{,}
\end{equation}
i.e. $e^{*}$ is non-singular uncorrelated and homoscedastic normal distributed. The residuals $e^{*}$ are called Theil's BLUS residuals.
In the short paper On the Testing of Regression Disturbances for Normality you find a comparison of OLS and BLUS residuals. In the tested Monte Carlo setting the OLS residuals are superior to BLUS residuals. But this should give you some starting point. 
