Expected value of Rayleigh quotient, non-centered Gaussian vector Let $X \sim \mathcal{N}\left(\mu, \Sigma \right)$, and let $A$ be a symmetric matrix. My understanding is that the Rayleigh quotient of vector $X$ is given by:
$$R=\frac{X^T A X}{X^T X}$$
I've been trying unsuccessfully to find an expression for the expected value of R. This seems to be a very well studied and relatively simple random variable, so I'd expect there to be some expression for $E(R)$.
I have seen these related questions(1, 2), but these ask about the generalized Rayleigh quotient, which is more complex.
I have also looked into resources discussing the distributions of quadratic forms. From his PhD thesis  beggining of section 3.2.2  (page 72), I understand that that $E(R) = \frac{E(X^T A X)}{E(X^T X)}$ (i.e. making the denominator $X^T I X$, since $I$ is idempotent as required), which seems to be off, and is in disagreement with simulations I've run. It is likely that I'm misunderstanding their claims though.
I also looked into Mathai & Provost "Quadratic forms in random variables", page 144 onwards (Ratios of Quadratic Forms), but couldn't find an answer. I may have missed the answer though, since most of the book goes over my head.
In case it's useful, this current question comes as a development of this previous question I asked here. It seems that a way to solve my previous problem is by estimating the expected value of the Rayleigh quotient.
Thanks

I found the answer for the case where the covariance of $X$ above is the identity ($\Sigma=I$). This was enough to solve my problem of interest, so I'm posting here and voting provisionally as the correct answer, until someone comes up with the more general case (which is related to the generalized Rayleigh quotient).
I found the answer to my question in this paper. Taking Equation 2.8 in that paper and substituting q=1* gives the following expectation:
$$\mathbb{E}\left[\frac{X^T A X}{X^T X}\right] = 
\frac{\operatorname{tr}(A)}{n} {}_1F_1\left(1; \frac{n}{2}+1; \frac{-\Vert\mu\Vert^2}{2}\right) +
\frac{1}{n+2} {}_1F_1\left(1; \frac{n}{2}+2; \frac{-\Vert\mu\Vert^2}{2}\right) \mu^T A \mu $$
for $X \sim \mathcal{N}(\mu, I)$. Here, $ {}_1F_1(a;b;c)$ is the confluent hypergeometric function, and $n$ is the number of dimensions in $X$.
I tested this expression in simulations and it seems to work perfectly. It is straightforward to extend it to a covariance that is $I\sigma^2$, where $\sigma^2$ is a scalar.
As a computational note, in case it's of use to anyone, the implementation of the confluent hypergeometric function in scipy was blowing up frequently for me. The package mpmath that allows for arbitrary precision computation has an implementation that has worked perfectly so far.
If you find this useful, consider upvoting the answer and/or question.
* The equation in the paper is more general, giving the expectation for $\mathbb{E}\left[\frac{X^T A X}{(X^T X)^q}\right]$
