# Distribution of the Rayleigh quotient

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.

• Do you know anything else about the relationship between $A$ and $B$ or properties beyond simply positive definiteness? Under a certain interpretation of "circular distribution" (i.e., invariant under orthogonal transformations), we can assume wlog that either $A$ or $B$ is diagonal. No assumption of positive definiteness of either is even necessary for that. Commented Jun 27, 2014 at 17:15
• A and B are correlation matrices. They are quite similar, but not identical. Commented Jun 27, 2014 at 19:01
• Possibly my choice of “circular distribution” was not ideal. What I mean is an elliptical distribution where the random variables w_i are independent – for example the standard normal distribution. Commented Jun 27, 2014 at 19:07

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.

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$$:

$$$$\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)$$$$

where: $$$$\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}$$$$

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