After some work I found the solution to the problem. The solution is the following:
Eq (1).
\begin{equation}
\mathbb{E}\left[\frac{XX^T}{X^TX}\right] =
\frac{1}{n}{}_{1}F_1\left(1; \frac{n}{2}+1; \frac{-||\mu||^2}{2}\right)I + \frac{1}{n+2} {}_{1}F_1\left(1; \frac{n}{2}+2; \frac{-||\mu||^2}{2}\right)\mu\mu^T
\end{equation}
for $X\sim\mathcal{N}(\mu, I)$, where $I$ in the formula is the identity matrix, $n$ is the dimensionality of $X$, and ${}_{1}F_1(a;b;c)$ is the confluent hypergeometric function. So, the solution is a weighted sum of the identity matrix and the outer product of the variable mean, with the weights given by the hypergeometric confluent function and dimension of $X$.
The derivation is as follows:
A well studied, related problem, is to find the expectation of ratios of quadratic forms $\mathbb{E}\left[\frac{X^TAX}{X^TX}\right]$, where $X\sim\mathcal{N}(\mu,I)$ and $A$ is a symmetric matrix. In this post, I had asked for an analytic formula to find this last expectation, and I then posted an answer I found. The formula is:
Eq (2).
\begin{equation}
\mathbb{E}\left[\frac{X^T A X}{X^T X}\right] =
\frac{tr(A)}{n} {}_1F_1\left(1; \frac{n}{2}+1; \frac{-||\mu||^2}{2}\right) +
\frac{1}{n+2} {}_1F_1\left(1; \frac{n}{2}+2; \frac{-||\mu||^2}{2}\right) \mu^T A \mu
\end{equation}
Now, we can use this formula to estimate each individual element $i,j$ of the original expectation by finding a matrix $A$ such that:
$$\mathbb{E}\left[\frac{XX^T}{X^TX}\right]_{i,j} =
\mathbb{E}\left[\frac{X_iX_j}{X^TX}\right] =
\mathbb{E}\left[\frac{X^T A X}{X^T X}\right]$$
It is easy to see that this relation is satisfied for each element $i,j$ by the following matrices:
- When $i=j$, $A_{i,i}=1$ and 0 elsewhere. In this case $tr(A)=1$ and $\mu^T A \mu = \mu_i^2$
- When $i \neq j$, $A_{i,j}=A_{j,i}=1/2$ and 0 elsewhere. In this case $tr(A)=0$ and $\mu^T A \mu = \mu_i\mu_j$
Using these matrices in Eq(2), we get, for $i=j$:
$$\mathbb{E}\left[\frac{XX^T}{X^T X}\right]_{i,i} =
\frac{1}{n} {}_1F_1\left(1; \frac{n}{2}+1; \frac{-||\mu||^2}{2}\right) +
\frac{1}{n+2} {}_1F_1\left(1; \frac{n}{2}+2; \frac{-||\mu||^2}{2}\right) \mu_i^2$$
and for $i\neq j$:
$$\mathbb{E}\left[\frac{XX^T}{X^T X}\right]_{i,j} =
\frac{1}{n+2} {}_1F_1\left(1; \frac{n}{2}+2; \frac{-||\mu||^2}{2}\right) \mu_i\mu_j$$
With a little of linear algebra, we notice that these two expressions equal the first formula in this answer.
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, please consider upvoting the answer and/or question.
