Confusion with matrix algebra when deriving GLS Consider the following linear model
$$y = X\beta + u$$
where $Var(u) = \mathbb{E}(uu') = \sigma^2\Omega$. I watched video on GLS derivation which proceeds as follows: since OLS is BLUE in case of homoskedasticity, we require such linear transformation (with transformation matrix $A$) of the regression equation that will bring us to this case, i.e. $Ay = AX\beta + Au$, where $Au = \sigma^2I$ and then estimate it with OLS. Since $u = y - X\beta$ we have the following:
\begin{equation}
\begin{aligned}
Var(Au) &= \sigma^2I \\
\Leftrightarrow \mathbb{E}(Auu'A') &= \sigma^2I \\
\Leftrightarrow A\mathbb{E}(uu')A' &= \sigma^2I \\
\Leftrightarrow \sigma^2AVA' &= \sigma^2I \\
\Leftrightarrow AVA' &= I  \\
\Leftrightarrow AV &= (A')^{-1} \\
\Leftrightarrow V &= A^{-1}(A')^{-1} \\
\text{(assuming $A$ is symmetric)} \; V &= A^{-2} \\
\text{(part of confusion)} \Leftrightarrow A &= V^{-1/2}
\end{aligned}
\end{equation}
I mean why the last step follows? First, $V$ is diagonal which simplifies the case. But we did not impose any structure on $A$ (except it is symmetric). $V^{-1/2}$ means take square root of each diagonal element. However, I am not convinced about the last step. In addition, would it work if I assume that $V$ is not diagonal? 
 A: It is not generally the case that the "matrix square root" just consists of the square roots of the elements on the main diagonal.
Consider a symmetric $V^{-1/2}$ via the Singular Value Decomposition (SVD), or eigenvalue decomposition, of $V$:
\begin{align*}
V&=PDP'\\
P'P&=PP'=I
\end{align*}
with $D$ diagonal. Then
\begin{eqnarray*}
V^{1/2}=PD^{1/2}P',
\end{eqnarray*}
and $V^{-1/2}$ is its inverse given by $PD^{-1/2}P'$.
Also, the "matrix square root" does not have to be symmetric (just as for conventional square roots, the "matrix square root" is not unique).
For example, for AR(1) errors $\epsilon_i=\phi\epsilon_{i-1}+\xi_i$, the variance-covariance matrix can be shown to be
$$
V=\sigma^2_{\xi}\frac{1}{1-\phi^2}\left(%
\begin{array}{cccc}
1 & \phi & \cdots & \phi^{n-1} \\
\phi & \ddots & & \phi^{n-2} \\
\vdots & \ddots & \ddots & \vdots \\
\phi^{n-1} & \cdots & \phi & 1 \\
\end{array}%
\right)
$$
One may verify that
$$
V(\phi)^{-1}=\frac{1}{\sigma^2_{\xi}}\left(%
\begin{array}{cccccc}
1 & -\phi & 0 & \cdots & \cdots & 0 \\
-\phi & 1+\phi^2 & -\phi & \ddots & & \vdots \\
0 & -\phi & \ddots & \ddots & \ddots & \vdots \\
\vdots & \ddots & \ddots & \ddots & \ddots & 0 \\
\vdots & & \ddots & \ddots & 1+\phi^2 & -\phi \\
0 & \cdots & \cdots & 0 & -\phi & 1 \\
\end{array}%
\right)
$$
One can also verify that the condition $V^{-1}=A'A$ is satisfied by the "Prais-Winsten transformation"
$$
A(\phi)=\frac{1}{\sigma_{\xi}}\left(%
\begin{array}{ccccc}
\sqrt{1-\phi^2} & 0 & \cdots & \cdots &0 \\
\hdashline -\phi & 1 & \ddots & & \vdots \\
0 & -\phi & 1 & \ddots & \vdots \\
\vdots & \ddots & \ddots & \ddots & 0 \\
0 & \cdots & 0 & -\phi & 1 \\
\end{array}%
\right)
$$ 
