Distribution of the Rayleigh quotient

Question

For a research project I need to find the expected value of the generalized Rayleigh quotient: $$E\,[w^T A w \ / \ w^T B w].$$ Here A and B are positive definite deterministic p x p covariance matrices, and w follows a multivariate distribution with circular altitude lines (say, multivariate standard normal). The dimension p is greater than 100.

This problem is easy to solve using simulation; however, I was wondering if somebody might know how this problem could be solved (or approximated) analytically. My first idea was that possibly by the Lindeberg or Lyapunov central limit theorem both the numerator and the denominator are approximately normal distributed, which gives us a ratio of two (correlated) normal random variables, but simulation shows that that is not the case.

Answer 1

In case of the normal distribution, a solution can be found in Mathai and Provost, Quadratic forms in random variables (1992). The inverse and product moments of such quadratic forms are derived there from the moment generating function.

Quadratic forms in elliptic distributions and their moments are treated in Mathai, Provost and Hayakawa, Bilinear forms and zonal polynomials (1995), but not to the same extend as in the normal case. As elliptical distributions are usually defined in terms of their characteristic function $e^{it\mu}\xi(t'\Sigma t)$, this function $\xi$ will appear in the solution if one chooses the mgf-approach. Yet, it has never been calculated, afaik.

Answer 2

There is a nice approximation described in the paper "Computing moments of ratios of quadratic forms in normal variables" (the approximation predates this paper though). It uses a second-order Taylor expansion that leads to a simple formula that is a good approximation in many cases (this approximation is used in this other answer of mine, see the comments of the original poster).

Let's write $N = w^T A w$ and $D = w^T B w$. Then $\mathbb{E}\left(\frac{w^T A w}{w^T B w}\right)$ can be approximated with the following expression of the moments of $N$ and $D$:

\begin{equation} \mathbb{E}\left(\frac{N}{D}\right) \approx \frac{\mu_N}{\mu_D}\left( 1 - \frac{Cov(N,D)}{\mu_N \mu_D} + \frac{Var(D)}{\mu_D^2} \right) \end{equation}

where: \begin{equation} \begin{split} & \mu_N = tr(A\Sigma) + \mu_{w}^T A \mu_{w} \\ & \mu_D = tr(B\Sigma) + \mu_w^T B \mu_w \\ & Var(D) = 2tr([B \Sigma]^2) + 4 \mu_w^T B \Sigma B \mu_w \\ & Cov(N,D) = 2tr(B \Sigma A \Sigma) + 4 \mu_w^T B \Sigma A \mu_w \end{split} \end{equation}

and $\mu_w$ and $\Sigma$ are the mean and covariance of normal vector $w$. That is, $w\sim \mathcal{N}(\mu_w, \Sigma)$.

If this answers your question, please consider upvoting.