Let $x$ and $y$ by unitvectors $\in\mathbb R^d$, such that $\|x\|_2=\|y\|_2=1$ and $\langle x,y\rangle=\rho\in[-1,1]$.
Let $G\in\mathbb R^{k\times d}$ be a matrix with independent Gaussian entries.
I'm interested in the mean and variance of $$\rho'=\frac{\langle Gx, Gy\rangle}{\|Gx\|_2\|G y\|_2}.$$
Another way to state this is that $Gx,Gy\in\mathbb R^k$ are Gaussian vectors with covariance matrix $\big(\begin{smallmatrix}1&\rho\\\rho&1\end{smallmatrix}\big)$. Thus $E[\langle Gx, Gy\rangle]=k \rho$ and $E[\|Gx\|_2^2]=E[\|Gy\|_2^2]=k$. From this we would expect $E[\rho']\approx\rho$.
I can make this more precise using the Gaussian Johnson Lindenstrauss lemma, which says that if $k=\varepsilon^{-2}\log1/\delta$, then $\langle Gx, Gy\rangle=\rho(1\pm\varepsilon)\|x\|_2\|y\|_2$ with probability at least $1-\delta$. We can get similar bounds for $\|Gx\|_2^2=\|x\|_2^2(1\pm\varepsilon)$ and for $y$. By Cauchy Schwartz we always have $\rho'\in[-1,1]$, so we get roughly
$$ E[\rho'] = \rho\,(1\pm O(1/\sqrt k)). $$
However, if I series expand $\rho'$ at $\rho=0$ I get
$$ \begin{align} E[\rho'] &= E\left[\frac{g}{\sqrt{\chi^2_{k-1} + g^2}} + \frac{\chi_{k} \chi^2_{k-1} }{(\chi^2_{k-1} +g^2)^{3/2}}\rho +O(\rho^2)\right] \\&= \frac{\Gamma(k/2+1/2)^2}{\Gamma(k/2+1)\Gamma(k/2)}\rho +O(\rho^2) \\&= \rho(1 + 1/(2k) + O(1/k^2)) +O(\rho^2), \end{align} $$ where $\chi_k^2$ and $\chi_{k-1}^2$ are independent Chi-Square distributed random variables, and $g$ is a standard Gaussian.
At the same time it seems $f(\rho)=E[\rho']$ for $\rho\ge 0$ is a convex function, that $f(\rho)\le\rho$, and that $f(1)=1$. Hence we must have $f(\rho)\ge \rho(1-1/(2k))$ and $$E[\rho'] = \rho(1+O(1/k)).$$
This is a quadratic improvement over the previous bound. I would love to have an intuition and a more rigorous proof. I also wonder what the variance of $\rho'$ is.