Hello! I have a question about how to calculate the expectation of a quadratic form as follows, where $X$ is a random variable that uniformly distributed on the unit sphere:
$$
E\_X[(\mathbf{x}^\top A\mathbf{x})^2]
=\int\_{\mathbf{x}\in S}{p(\mathbf{x})(\mathbf{x}^\top A\mathbf{x})^2dS(\mathbf{x})}
=\int\_{\mathbf{x}\in S}{\frac{1}{4\pi}(\mathbf{x}^\top A\mathbf{x})^2dS(\mathbf{x})}
$$
where $S=\{\mathbf{x}\in\mathcal{R}^N|\mathbf{x}^\top\mathbf{x}=1\}$.


If $\mathbf{x}$ is Gaussian, there are some conclusions about the expectation of the quadratic forms. However, I find it difficult to deal with when the variable is distributed on a sphere. When the dimension is 2 or 3, this problem can be solved by representing it with polar coordinate. However, when the dimension is high, such a representation will be rather redundant, how can I calculate this integral then?
Please give some help for this problem if you have any idea. Thank you very much!



How to calculate this expectation where the random variable is restricted on a sphere?
Assuming $A$ is diagonalizable in an orthonormal basis with a diagonal of $a\_k$'s, one asks for the expectation $C(A)$ of the random variable
$$
(X^TAX)^2=\left(\sum\limits\_ka\_kX\_k^2\right)^2=\sum\limits\_ka\_k^2X\_k^4+\sum\limits\_{k\ne \ell}a\_ka\_\ell X\_k^2X\_\ell^2,
$$
where the vector $X=(X\_k)\_{1\leqslant k\leqslant d}$ is uniformly distributed on the Euclidean unit sphere. Hence,
$$
C(A)=c\_d\cdot\sum\_ka\_k^2+b\_d\cdot\sum\_{k\neq\ell}a\_ka\_\ell=(c\_d-b\_d)\cdot\sum\_ka\_k^2+b\_d\cdot\left(\sum\_ka\_k\right)^2,
$$
with
$$
c\_d=\mathrm E(X\_1^4),\qquad b\_d=\mathrm E(X\_1^2X\_2^2).
$$
Since the $X\_k^2$'s are identically distributed and sum to $\|X\|^2=1$, $\mathrm E(X\_1^2)=1/d$ and
$$
dc\_d+d(d-1)b\_d=1.
$$
Furthermore, a direct computation **(1)** yields $b\_d=1/(d(d+2))$ hence $c\_d-b\_d=2/(d(d+2))$ and
$$
C(A)=\frac1{d(d+2)}\cdot\left(2\sum\limits\_ka\_k^2+\left(\sum\limits\_ka\_k\right)^2\right).
$$
Finally,
$$
\color{red}{\mathrm E((X^TAX)^2)=\frac{2\text{tr}(A^2)+\left(\text{tr}(A)\right)^2}{d(d+2)}}.
$$
**Note:**


**(1)** One can write $X$ as $X=Z/\|Z\|$ where the vector $Z=(Z\_k)\_{1\leqslant k\leqslant d}$ is i.i.d. standard gaussian. Thus,
$$
X\_1^2X\_2^2=Z\_1^2Z\_2^2\|Z\|^{-4}=Z\_1^2Z\_2^2\,\int\_0^{+\infty}\mathrm e^{-t\|Z\|^2}\cdot t\mathrm dt,
$$
hence
$$
X\_1^2X\_2^2=\int\_0^{+\infty}Z\_1^2\mathrm e^{-tZ\_1^2}\cdot Z\_2^2\mathrm e^{-tZ\_2^2}\cdot \mathrm e^{-t(Z\_3^2+\cdots+Z\_d^2)}\cdot t\mathrm dt.
$$
Integrating this,
$$
b\_d=\mathrm E(X\_1^2X\_2^2)=\int\_0^{+\infty}u'(t)^2\cdot u(t)^{d-2}\cdot t\mathrm dt\quad\text{with}\quad u(t)=\mathrm E(\mathrm e^{-tZ\_1^2}).
$$
A standard computation yields $u(t)=(1+2t)^{-1/2}$, hence $u'(t)=-(1+2t)^{-3/2}$ and
$$
b\_d=\int\_0^{+\infty}(1+2t)^{-d/2-2}\cdot t\mathrm dt=\left.\frac1{2d(1+2t)^{d/2}}-\frac1{2(d+2)(1+2t)^{d/2+1}}\right|\_0^{+\infty}=\frac1{d(d+2)}.
$$


