Solution (source code)

= Solution

By the <spectral theorem>, write $V=U\operatorname{diag}(v_1,\ldots,v_p)U^T$, where $U$ is <orthogonal> and each $v_j>0$. Choose the <whitening transformation>
$$
\boxed{L=V^{-1/2}=U\operatorname{diag}(v_1^{-1/2},\ldots,v_p^{-1/2})U^T.}
$$
A <linear image of a multivariate normal vector> is <multivariate normal>, and $LVL^T=I$. Thus $Z_i=L(Y_i-\mu)$ are independent $N_p(0,I)$ <random vectors>. Their whole joint distribution is independent of $\mu,V$.

Define $\bar Z=n^{-1}\sum_iZ_i$ and $S_Z=n^{-1}\sum_i(Z_i-\bar Z)(Z_i-\bar Z)^T$. These satisfy
$$
\bar Z=L(\bar Y-\mu),\qquad S_Z=LSL^T,
\qquad S_Z^{-1}=L^{-T}S^{-1}L^{-1}.
$$
Therefore
$$
\boxed{(\bar Y-\mu)^TS^{-1}(\bar Y-\mu)
=\bar Z^TS_Z^{-1}\bar Z.}
$$
The right side is a function of independent standard <multivariate normal> observations only, so its distribution depends on $n,p$ but not on $\mu,V$. This is the same <affine invariance of Hotelling's statistic> that permits multivariate mean tests with unknown <covariance>. The proof concerns the jointly transformed mean and <covariance>; whitening the mean alone would not establish the conclusion. Inverses exist almost surely under $n>p$. If $V$ is singular, no square <matrix> can whiten it to $I$, so nonsingularity of the population <covariance> is necessary as well.