variance of the square of the bias on linear regression

Question

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

Answer 1

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:

1. 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.

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

3. 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).$$

Answer 2

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). $$