# Variance of quadratic form for multivariate normal distribution

This is a homework problem I’m trying to solve but I can’t seem to solve Q1b without using the theorem. I am also given the fact that

$$E(y’Ay)=tr(A\Sigma)+\mu’A\mu$$

I’ve tried using the trace-expectation trick but to no avail. Assuming $$y\sim N(0,I)$$, \begin{align*} var(y’Ay) &=E[(y’Ay-tr(A))^2]\\ &=E(y’Ayy’Ay)-tr(A)^2\\ &=E(tr[y’Ayy’Ay])-tr(A)^2\\ &=E(tr[Ayy’Ayy’])-tr(A)^2\\ &=tr(E[Ayy’Ayy’])-tr(A)^2\\ &=tr(AE[yy’Ayy’])-tr(A)^2 \end{align*} Then I’m stuck.

Suppose $$A=((a_{ij}))$$ and $$y=(y_1,y_2,\ldots,y_n)'\sim N(0,I_n)$$, so that $$y_i$$s are i.i.d standard normal.

You already have $$\operatorname E[y'Ay]=\operatorname{tr}(A)$$ by the result you quoted (discussed here).

For the variance, you can simply use $$\operatorname{Var}(y'Ay)=\operatorname E[(y'Ay)^2]-(\operatorname E[y'Ay])^2\tag{1}$$

To compute the first expectation, we write the quadratic form as $$y'Ay=\sum_{i,j}a_{ij}y_iy_j$$ to get $$(y'Ay)^2=\sum_{i,j,k,l}a_{ij}a_{kl}y_iy_jy_ky_l \tag{2}$$

Now observe that $$\operatorname E[y_iy_jy_ky_l]=\begin{cases}3&,\text{ if }i=j=k=l \\ 1&,\text{ if }i=j,k=l;i=k,j=l;i=l;j=k \\ 0&,\text{ otherwise }\end{cases}$$

Therefore taking expectation on both sides of $$(2)$$,

$$\operatorname E[(y'Ay)^2]=3\sum_i a_{ii}^2+\sum_i\left(\sum_{k\ne i}a_{ii}a_{kk}+\sum_{j\ne i}a_{ij}^2+\sum_{j\ne i}a_{ij}a_{ji}\right)\tag{3}$$

Keeping in mind that $$A$$ is symmetric, you have $$\operatorname{tr}(A^2)=\sum_{i,j}a_{ij}^2$$.

It is now straightforward to see that $$(3)$$ reduces to $$\operatorname E[(y'Ay)^2]=(\operatorname{tr}(A))^2+2\operatorname{tr}(A^2)$$

From $$(1)$$ you get the desired result $$\boxed{\operatorname{Var}(y'Ay)=2\operatorname{tr}(A^2)}$$

You can now try to generalize this to the case $$y\sim N(0,\Sigma)$$ and hence to $$y\sim N(\mu,\Sigma)$$.

Reference:

• Linear Regression Analysis by Seber and Lee.
• Since $A$ is symmetric, the easier way to find the variance directly for the $N(0,I_n)$ case is to use spectral decomposition to write $A=P'\Lambda P$ for some orthogonal $P$ and $\Lambda=\operatorname{diag}(\lambda_1,\ldots,\lambda_n)$ where $\lambda_i$ are eigenvalues of $A$. Then $y'Ay=x'\Lambda x=\sum_{i=1}^n \lambda_i x_i^2$ with $x=Py\sim N(0,I_n)$. Jul 14 '20 at 16:12