product of quadratic forms of random vectors uniform on the sphere Let $g = (g_1, ..., g_n)$ be a random vector distributed uniformly on the sphere $\{ x \in \mathbb{R}^n : \| x \|_2 = 1 \}$. Let $A, B$ be two symmetric $n \times n$ matrices. 
I am interested in a simple formula for: $$\mathbb{E}[ g^T A g g^T B g ]$$
In particular, I know that if $g$ is a standard normal $N(0, I)$, then 
$$
  \mathbb{E}[ g^T A g g^T B g ] = 2 \mathrm{Tr}(AB) + \mathrm{Tr}(A)\mathrm{Tr}(B) \:.
$$
Does something similar hold in the uniform on a sphere case?
 A: It turns out the answer is:
$$
\mathbb{E}[g^T A g g^T B g] = \frac{1}{n(n+2)} (2 \mathrm{Tr}(AB) + \mathrm{Tr}(A) \mathrm{Tr}(B)) \:.
$$
This comes from the Theorem in Section 3 of 
D. P. Wiens, "On Moments of Quadratic Forms in Non-Spherically Distributed Variables" (https://pdfs.semanticscholar.org/f9e3/9424ff32aec189b441dd189b6f819764de6b.pdf). It is straightforward to check that the conditions (1.1) hold for the uniform distribution over the sphere. 
To apply the result, you need to compute $\mathbb{E}[g_1^2 g_2^2]$. It turns out this is equal to $\frac{1}{n(n+2)}$. Here is a simple way to compute this by symmetry. Observe that:
$$
g_1^2 g_2^2 = (1 - g_2^2 - ... - g_n^2) g_2^2 = g_2^2 - g_2^4 - g_3^2 g_2^2 - ... - g_n^2 g_2^2 \:.
$$
Let $\mu = \mathbb{E}[g_1^2 g_2^2]$. By symmetry, $\mu = \mathbb{E}[g_i^2 g_j^2]$ for any $i \neq j$. Hence taking expectations,
$$
 \mu = \mathbb{E}[g_2^2] - \mathbb{E}[g_2^4] - (n-2) \mu \:.
$$
Rearrange to obtain:
$$
\mu = \frac{1}{n-1}\left( \mathbb{E}[g_2^2] - \mathbb{E}[g_2^4] \right) \:.
$$
We can compute $\mathbb{E}[g_2^2] = 1/n$ by a similar symmetry argument. Furthermore, since $g_2^2 \stackrel{d}{=} \mathrm{Beta}(1/2, (n-1)/2)$, we can compute $\mathbb{E}[g_2^4]$ by looking up the formula of the second moment of a Beta distribution, which gives us $\mathbb{E}[g_2^4] = \frac{3}{n(n+2)}$. Hence,
$$
\mu = \frac{1}{n-1}\left( \frac{1}{n} - \frac{3}{n(n+2)} \right) = \frac{1}{n(n+2)} \:.
$$
A: Here is another approach to reduce the problem to the Gaussian case. Let $g \sim N(0,I)$ and let $\theta = g / \|g\|$ and $r = \|g\|$ so that $g = r\theta$. It is not hard to verify that $\theta$ is uniformly distributed on the sphere and is independent of $r$. Then, we have, by independence,
\begin{align*}
  \mathbb E (g^T A g g^T B g) = \mathbb E (\theta^T A \theta \theta^T B \theta) \cdot \mathbb E (r^4)
\end{align*}
and
\begin{align*}
  \mathbb E (r^4) &= \mathbb E \Big(\sum_i g_i^2\Big)^2 \\
    &=\sum_i \mathbb E  g_i^4 +  \sum_{i \neq j} \mathbb E g_i^2 \mathbb Eg_j^2   \\
&= 3 n + n(n-1) = n(n+2).
\end{align*}
