Let $S^{d-1}$ denote the $d$-dimensional sphere, and $\sigma_{d-1}$ denote the corresponding spherical measure. I am wondering how to go about solving the following integral for $c>0$
$$ \int_{S^{d-1}} \|M u\|_2^{2c} d \sigma_{d-1}(u), $$
where $M$ is a positive definite and symmetric matrix. The obvious first step is to substitute $v=Mu$ but I am not quite sure how to work with the spherical measure, even for the simpler setting of $c=1$. In the setting where $M=I$ then $\|Mu\|= \|u\|=1$ whenever $u \in S^{d-1}$, so here we have $$ \int_{S^{d-1}} \|M u\|_2^{2c} d \sigma_{d-1}(u) = \int_{S^{d-1}} \|u\|_2^{2c} d \sigma_{d-1}(u) = \int_{S^{d-1}} d \sigma_{d-1}(u) = S(d) $$
where $S(d)$ is the surface area of the $d$-dimensional sphere.
Any hints or references to relevant material for dealing with the more general setting would be appreciated!
Actually $v=Mu$ is the wrong way to go, as you said you don't know how to work with the spherical measure (it doesn't even preserve the sphere).
The trick is that since $M$ is symmetric and positive definite, we have an orthogonal matrix $O$ (that is, $O^T O = O O^T=I$) such that $$ M = O^T \Lambda O, $$ where $\Lambda = \operatorname{diag}(\lambda_1,\dots,\lambda_d)$ is a diagonal matrix with positive entries.
Then $$ \|Mu\|_2^2=u^TM^TMu=u^TO^T \Lambda OO^T \Lambda Ou=(Ou)^T\Lambda^2(Ou)=\|\Lambda(Ou)\|_2^2. $$ Now let $v=Ou$, and since $O$ is orthogonal, it is a diffeomorphism of the sphere, and $d\sigma_{d-1}(u)=d\sigma_{d-1}(v)$. We then have $$ \int_{S^{d-1}} \|Mu\|_2^{2c}d\sigma_{d-1}(u)=\int_{S^{d-1}} \|\Lambda(Ou)\|_2^{2c}d\sigma_{d-1}(u)=\int_{S^{d-1}} \|\Lambda v\|_2^{2c}d\sigma_{d-1}(v). $$
Now $$ \|\Lambda v\|_2^2 = \lambda_1^2 v_1^2 + \dots + \lambda_d^2 v_d^2. $$ This integral is easier than before. I am not too sure how to go for a real $c>0$, but if $c$ is a positive integer, we can expand $(\|\lambda v\|_2^2)^c$ out, and reduce it to a sum of the following kind of integrals.
We actually have a concrete formula $$ \int_{S^{d-1}} v_1^{2a_1}\dots v_d^{2a_d}\, d\sigma_{d-1}(v) = \frac{2\Gamma(a_1+\frac{1}{2})\cdots \Gamma(a_d+\frac{1}{2})}{\Gamma(a_1+\cdots+a_d+\frac{d}{2})}, \quad a_i\in {\mathbb Z}_{\geq 0}. $$ This can be standard, but let me write out the details for the proof.
Consider the spherical coordinates on ${\mathbb R}^d$ with $\rho>0$, $(v_1,\dots,v_d)\in S^{d-1}$. Then $(x_1,\dots,x_d)\in {\mathbb R}^d$ with $x_i=\rho v_i$. We have \begin{align} &\quad\ \big(\int_{S^{d-1}} v_1^{2a_1}\dots v_d^{2a_d}\, d\sigma_{d-1}(v)\big)\big(\int_0^\infty e^{-\rho^2} \rho^{2(a_1+\dots+a_d)+d-1}\,d\rho\big)\\ &=\int_0^\infty e^{-\rho^2} \int_{S^{d-1}} (\rho v_1)^{2a_1}\dots (\rho v_d)^{2a_d} d\sigma_{d-1}(v)\, \rho^{d-1}d\rho\\ &=\int_{{\mathbb R}^d} e^{-(x_1^2+\dots+x_d^2)}x_1^{2a_1}\dots x_d^{2a_d}\,dx\\ &=\big(\int_{\mathbb R} e^{-x_1^2}x_1^{2a_1}\,dx_1\big)\dots \big(\int_{\mathbb R} e^{-x_d^2}x_d^{2a_d}\,dx_d\big). \end{align} Now simple substitutions show that \begin{align} \int_0^\infty e^{-\rho^2} \rho^{2(a_1+\dots+a_d)+d-1}\,d\rho&= \frac{1}{2}\Gamma(a_1+\dots+a_d+\frac{d}{2}),\\ \int_{\mathbb R} e^{-x_i^2}x_i^{2a_i}\,dx_i &= \Gamma(a_i+\frac{1}{2}). \end{align} Putting these together, we get the desired formula.