variance of the square of the bias on linear regression Basic setting
let the linear model be:
$$
\mathbf{y}=\mathbf{X\beta}+\epsilon
$$
where $\epsilon \sim N(0,\sigma^2\mathbf{I}_n)$
$n$ is the number of samples
$p$ is the number of attributes.
$\mathbf{y}\in\mathbb{R}^{n \times 1}$, is known.
$\mathbf{X}\in\mathbb{R}^{n \times p}$, is known.
$\mathbf{\beta}\in\mathbb{R}^{p \times 1}$, is unknown.
we estimate $\beta$ by minimizing the least squares,and we have:
$$
\hat \beta = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}\mathbf{y}
$$
question
let
$$
L_1^2 = (\hat \beta - \beta)^T(\hat \beta - \beta)
$$
show that
$$
Var(L_1^2)=2 \sigma^4 \text{Trace}((X^TX)^{-2})
$$ 
What I have known:
$$
Var(\hat \beta)= \sigma^2 (X^TX)^{-1}
$$
$$
E(L_1^2)=\sigma^2 \text{Trace}((X^TX)^{-1})
$$ 
I meet this question when I was reading Ridge regression: Biased estimation for nonorthogonal problems Hoerl, Arthur E;Kennard, Robert W Technometrics; Feb 2000; 42, 1; ProQuest pg. 80
 A: You did the hard part already.
Let's simplify notation a little.  Notice
$$\hat\beta - \beta = (X^\prime X)^{-1}X^\prime (X\beta y + \epsilon) - \beta = (X^\prime X)^{-1}X^\prime \epsilon.$$
Therefore we may write $L_1^2$ as
$$(\hat\beta-\beta)^\prime(\hat\beta-\beta)= \epsilon^\prime X^\prime (X^\prime X)^{-2} X \epsilon = \epsilon^\prime A \epsilon = \sum_{i,j} \epsilon_i\, a_{ij}\, \epsilon_j.$$
Note that $A$ is symmetric: $a_{ij} = a_{ji}$ for all indexes $i$ and $j.$  Moreover,
$$\operatorname{Tr}(A) = \operatorname{Tr}\left(X^\prime (X^\prime X)^{-2} X\right)=\operatorname{Tr}\left(X^\prime X(X^\prime X)^{-2} \right)=\operatorname{Tr}\left((X^\prime X)^{-1}\right)$$
and similarly 
$$\operatorname{Tr}(A^2) = \operatorname{Tr}\left((X^\prime X)^{-2}\right).$$
Choose units of measurement for the $y_i$ that make $\sigma^2=1$ so we don't have to track it: we know this will introduce a factor of $\sigma^4$ at the end.
The only fact about Normal variates we will need is that when the $\epsilon_i$ are independent standard Normal variables,
$$E[\epsilon_i\epsilon_j\epsilon_k\epsilon_l] = \delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{kj}$$
where $\delta_{ij} = 1$ when $i=j$ and $0$ otherwise is the Kronecker delta.  This scarcely needs proof, because a little reflection on its structure shows it merely states the following:


*

*The expectation is zero unless the $\epsilon$'s can be paired up, because otherwise the symmetry of the standard Normal distribution shows the expectation equals its negative.

*When two of the $\epsilon$'s are equal, they introduce a factor of $1$ in the expectation (because they have unit variance).

*In the special case where all four of the $\epsilon$'s are equal, we obtain the Normal kurtosis, which is $3.$
To compute the variance, we need to find the expected square, which is accomplished by invoking the foregoing result and the linearity of expectation:
$$\eqalign{
E[((\hat\beta-\beta)^\prime(\hat\beta-\beta))^2] &= E\left[\sum_{i,j}\epsilon_i\, a_{ij}\, \epsilon_j\ \sum_{k,l}\epsilon_k\, a_{kl}\, \epsilon_l\right] \\
&= \sum_{i,j,k,l} a_{ij} a_{kl} \left(\delta_{ij}\delta_{kl} + \delta_{ik}\delta_{jl} + \delta_{il}\delta_{kj}\right) \\
&= \sum_{i,k} a_{ii}a_{kk} + \sum_{i,j}a_{ij}a_{ij} + \sum_{i,k}a_{ik} a_{ki} \\
&=\operatorname{Tr}(A)^2 + 2\operatorname{Tr}(A^2).
}$$
Subtracting off $(E[L_1^2])^2 = \operatorname{Tr}(A)^2$ yields the variance which--in terms of the original unit of measure $\sigma$--is
$$\operatorname{Var}(L_1^2) = 2\sigma^4\operatorname{Tr}(A^2) = 2\sigma^4\operatorname{Tr}\left((X^\prime X)^{-2}\right).$$ 
A: From the assumption that $\varepsilon \sim \mathcal N(0,\sigma^2 I)$ we have
$$
\hat\beta \sim \mathcal N(\beta, \sigma^2 (X^TX)^{-1})
$$
so
$$
\gamma := \hat\beta - \beta \sim \mathcal N(0, \sigma^2 (X^TX)^{-1}).
$$
We now are trying to get the variance of the Gaussian quadratic form $\gamma^T\gamma = L_1^2$. I'll do this by working out the moment generating function of $\gamma^T\gamma$.
$$
\text E(e^{t\gamma^T\gamma}) = \int e^{t\gamma^T\gamma} \frac{\vert X^TX\vert^{1/2}}{(2\pi\sigma^2)^{p/2}}\exp\left(-\frac 1{2\sigma^2}\gamma^TX^TX\gamma \right)\,\text d\gamma \\
= \frac{\vert X^TX\vert^{1/2}}{(2\pi\sigma^2)^{p/2}} \int \exp\left(-\frac 1{2\sigma^2}\gamma X^TX \gamma + t\gamma^T\gamma\right)\,\text d\gamma.
$$
Inside the exponential we have
$$
-\frac 1{2\sigma^2}\gamma X^TX \gamma + t\gamma^T\gamma = -\frac 1{2\sigma^2}\gamma^T \left[X^TX - 2\sigma^2 t I\right]\gamma
$$
and $X^TX$ being invertible means $\lambda_{\min}(X^TX) > 0$ and for $t > 0$ sufficiently small we'll have $\lambda_{\min}(X^TX - 2\sigma^2 t I) > 0$ too which means that there is a $\delta > 0$ such that $t \in (0, \delta) \implies X^TX - 2\sigma^2 t I$ is invertible. And for $t \leq 0$ this matrix is invertible too so this means $M(t)$ is finite on an interval containing $0$ which means it's safe to use moment generating functions for this.
Taking $t$ to be sufficiently small, we have
$$
M(t) := \text E(e^{t \gamma^T\gamma}) =  \frac{\vert X^TX\vert^{1/2}}{(2\pi\sigma^2)^{p/2}} \cdot (2\pi\sigma^2)^{p/2} |X^TX - 2\sigma^2 t I|^{-1/2} \\
= \frac{|X^TX|^{1/2}}{|X^TX - 2\sigma^2 t I|^{1/2}}.
$$
Now we can use some matrix calculus (all the results that you need for this are in the various tables in the wikipedia article on matrix calculus) to find
$$
M'(t) = -\frac 12 |X^TX|^{1/2}|X^TX - 2\sigma^2 t I|^{-3/2} \cdot |X^TX - 2\sigma^2 t I| \cdot \text{tr}\left((X^TX - 2\sigma^2 t I)^{-1} (-2 \sigma^2 I)\right)
$$
so
$$
M'(0) = -\frac 12 |X^TX|^{1/2}|X^TX|^{-3/2}|X^TX|\text{tr}\left((X^TX)^{-1} (-2 \sigma^2 I)\right) \\
= \sigma^2 \text{tr}((X^TX)^{-1})
$$
which confirms the first moment (although if this is all we were doing it'd have been far easier to just use the usual trick of switching traces and expectations -- there may be a similar trick for the variance but if there is I don't know it).
Now for the second moment, we can go through some more matrix calculus to get
$$
M^{\prime\prime}(t) = \sigma^2 |X^TX|^{1/2} \frac{\partial }{\partial t} \left[|X^TX - 2\sigma^2 t I|^{-1/2} \cdot \text{tr}\left((X^TX - 2\sigma^2 t I)^{-1} \right)\right] \\
= \sigma^2 |X^TX|^{1/2} \left[ -\frac 12 \cdot |X^TX - 2\sigma^2 t I|^{-3/2} \cdot |X^TX - 2\sigma^2 t I| \cdot \text{tr}\left((X^TX - 2\sigma^2 t I)^{-1} (-2 \sigma^2 I)\right)^2 \\
+ |X^TX - 2\sigma^2 t I|^{-1/2} \cdot 2\sigma^2 \text{tr}\left((X^TX - 2\sigma^2 t I)^{-2}\right)\right]
$$
so
$$
M^{\prime\prime}(0) = \sigma^2 |X^TX|^{1/2} \left[\sigma^2 |X^TX|^{-1/2}\text{tr}\left((X^TX)^{-1} )\right)^2 + 2\sigma^2 |X^TX|^{-1/2}\text{tr}\left((X^TX)^{-2}\right)\right] \\
= \sigma^4 \text{tr}\left((X^TX)^{-1} )\right)^2 + 2\sigma^4 \text{tr}\left((X^TX)^{-2}\right)
$$
so
$$
\text{Var}(L_1^2) = M''(0) - M'(0)^2 = 2\sigma^4 \text{tr}\left((X^TX)^{-2}\right).
$$
