How to calculate this expectation where the random variable is restricted on a sphere? 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!
 A: 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)}.
$$
A: See http://arxiv.org/abs/math/0305252 for computations of this sort.
