Stuck applying the law of the unconscious statistician I am trying to compute the expectation the matrix $g(x):=xx^\top / x^\top H x$ with $x$ a multivariate normal $\mathcal{N}(0,\Sigma)$ and $H$ fixed positive definite, but I am stuck when integrating by parts.
I thought to focus on the $ij$-th element, i.e., using the law of the unconscious statistician,
$$
\mathbb{E}_u\left[
\frac{x_i x_j}{\sum_{k,\ell} x_k x_l H_{k\ell}}
\right]
=\int_{-\infty}^\infty
\underbrace{
\frac{x_i x_j}{\sum_{k,\ell} x_k x_l H_{k\ell}}
}_{g(x)}
\cdot
\underbrace{
\frac{\exp\left(-\frac 1 2 x^\top\Sigma^{-1}x\right)}{
\sqrt{(2\pi)^d|\Sigma|}
}
}_{p(x)}
dx
$$
Integration by parts is, in general:
$$
\int_{-\infty}^\infty u(x)v(x)dx
=\left[u(x)\int v(x)dx\right]_{-\infty}^\infty
-\int_{-\infty}^\infty \left( du(x) \int v(x) dx \right)
$$
Now I must use $u:=p$ and $v:=g$, because the other way would result in $\int p(x)dx$ for which no elementary solution exists.
Next I tried to solve $\int g(x)dx$ in the simpler two-dimensional case with $i=1$ and $j=2$, i.e.,
$$
\int\frac{x_1 x_2}{x_1^2 h_{11}+x_2^2 h_{22}+2x_1x_2h_{12}}dx_1dx_2
$$
Plugging this into sympy gives a very long and ugly expression that I don't think you want to see :) and I don't think this path will lead anywhere.
How would you proceed to compute or approximate this expectation?
 A: Suppose $X \in \mathbb{R}^d$. Applying the spectral theorem, we may write $\Sigma^{1/2}\mathbf{H}\Sigma^{1/2} = \mathbf{U}^{\top}\mathbf{D}\mathbf{U}$ for an orthogonal matrix $\mathbf{U}$ and a diagonal matrix $\mathbf{D} = \operatorname{diag}(\lambda_1,\ldots,\lambda_d)$. Then the distribution of $X$ can be realized as
$$X \sim \mathbf{A} R \Theta,$$
where $\mathbf{A} = \Sigma^{1/2} \mathbf{U}^{\top}$, and $R \sim \chi(d)$ and $\Theta\sim\operatorname{Uniform}(\mathbb{S}^{d-1})$ independent random variables. So,
$$ \frac{XX^{\top}}{X^{\top}\mathbf{H}X} \sim \mathbf{A} \frac{\Theta\Theta^{\top}}{\Theta^{\top}\mathbf{D}\Theta} \mathbf{A}^{\top}. $$
Taking expectation to both sides,
$$ \mathbb{E}\left[\frac{XX^{\top}}{X^{\top}\mathbf{H}X}\right]
= \mathbf{A} \mathbf{G} \mathbf{A}^{\top}, \qquad \mathbf{G} = \mathbb{E}\left[ \frac{\Theta \Theta^{\top}}{\Theta^{\top}\mathbf{D}\Theta} \right]. $$
Now by writing the $i$th component of $\Theta$ as $\Theta_i$, the $(i,j)$-component of $\mathbf{G}$ is given by
$$ G_{ij}
= \mathbb{E}\left[ \frac{\Theta_i \Theta_j}{\sum_k \lambda_k \Theta_k^2} \right]. $$

*

*By the reflection symmetry, we have $G_{ij} = 0$ whenever $i \neq j$.


*When $d = 2$, we have $G_{11} = \frac{1}{\sqrt{\lambda_1}(\sqrt{\lambda_1} + \sqrt{\lambda_2})}$ and $G_{22} = \frac{1}{\sqrt{\lambda_2}(\sqrt{\lambda_1} + \sqrt{\lambda_2})}$.


*However, I guess that no such closed-form exact formula exists for high dimensions. Indeed, even for $d = 3$, the formula quickly goes out of hands.


*On the other hand, when $d \gg 1$ and $\lambda_i \ll 1 \ll \sum_{k=1}^{d} \lambda_k$, LLN suggests that $G_{ii} \approx \frac{1}{\sum_{k=1}^{d} \lambda_k} $.


*Also note that $\operatorname{Tr}(\mathbf{D}\mathbf{G}) = \sum_{i=1}^{n} \lambda_i G_{ii} = 1$.
A: Here is a variation on a prior answer.
First change variables so that (i) $X^T HX= Y^T Y$ and (ii)  $X^T\Sigma X = Y^T \Lambda Y$ where $\Lambda $ is diagonal and positive definite. This is possible by the usual diagonalization theory for a pair of  symmetric quadratic forms (e.g. spectral theorem.)
Then you can use odd-symmetry to see that in these new $Y$ variables, the expectations of non-diagonal terms such as $\int \frac{ y_i y_j}{||Y||^2} e^{-\sum_i \lambda_i y_1^2}$ are zero.
For the diagonal term integrals $M(\lambda_1, \lambda_2, \ldots, \lambda_n)=\int \frac{ y_i^2}{||Y||^2} e^{-\sum_i \lambda_i y_1^2}$ , I am out of ideas. Maybe try differentiating the unknown value of the integral $M(\lambda_1, \lambda_2, \ldots, \lambda_n)$ w.r.t. to the sum of all the various $\lambda_i$ to kill off denominator and obtain standard Gaussian moment integrals, to get a derivative identity that can be unraveled?
A: I hope that someone finds a mathematical approximation to answer this question, but for the present you can quickly simulate the result for a particular choice of multivariate Gaussian distribution and positive-semi-definite matrix $H$.
import numpy as np
from sklearn.datasets import make_spd_matrix
import matplotlib.pyplot as plt

np.random.seed(0)

n = 3
m = 1000
mu = np.random.normal(0,10,size=n)
Sigma = make_spd_matrix(n)
H = make_spd_matrix(n)


v = np.random.multivariate_normal(mu, Sigma, size=m)
instances = []
for i in range(m):
    for j in range(i,m):
        instances.append(v[i] @ v[j] / (v[i] @ H @ v[j]))

plt.hist(instances, bins=100)
plt.title(f'mean:{np.mean(instances)}\nStd:{np.std(instances)}')
plt.show()


