Average value of a bilinear map on a Euclidean sphere

I figured it out, here goes:

Let $(e_1, \dots, e_n)$ be basis of $V$ which is $g$-orthonormal and $B$-orthogonal (the existence of such a basis is precisely the spectral theorem). Let $\lambda_k = B(e_k, e_k)$ for $k \in \{1, \dots n\}$. For any vector $x = \sum_{k=1}^n x_k e_k$, the squared norm of $x$ is $g(x,x) = \sum_{k=1}^n {x_k}^2$ and the quadratic form is given by $B(x,x) = \sum_{k=1}^n \lambda_k {x_k}^2$. Therefore: \begin{equation} \int_{S_r} B(x, x) \sigma_r = \sum_{k=1}^n \lambda_k \int_{S_r} {x_k}^2 \sigma_r~. \end{equation} It remains to compute the integrals $I_k = \int_{S_r} {x_k}^2 \sigma_r$. This can be done swiftly using a little trick, starting with the observation that any two of these integrals are equal. Indeed, for any $k \neq l$, one can easily find an isometry $\varphi \in O(g)$ such that ${x_k}^2 \circ \varphi = {x_l}^2$ (namely, the orthogonal reflection through the line spanned by $e_k + e_l$). Since $\varphi$ preserves $\sigma_r$, the change of variables theorem ensures that $I_k = I_l$. Since all the integrals $I_k$ are equal, one can write $I_k = \frac{1}{n} \sum_{l=1}^n I_l$ for any $k$. That is $I_k = \frac{1}{n}\int_{S_r} \left(\sum_{k=1}^n {x_k}^2\right) \sigma_r$. However $\sum_{k=1}^n {x_k}^2 = g(x,x) = r^2$ for any $x \in S_r$. This yields $I_k = \frac{1}{n}\int_{S_r} r^2 \sigma_r = \frac{r^2}{n} \operatorname{vol}(S_r)$. Coming back to the initial integral, we obtain the desired result: \begin{equation} \int_{S_r} B(x, x) \sigma_r = \frac{r^2}{n} \operatorname{vol}(S_r) \sum_{k=1}^n \lambda_k = \frac{r^2}{n} \operatorname{vol}(S_r) \operatorname{tr}_g (B)~. \end{equation}