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} $$


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


Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.