# 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

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

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