Distribution of the norm of uniform random unit vector after linear transformation

Question

Suppose that $\mathbf{u}$ is a uniform unit vector. It is obtained as $\mathbf{u}=\frac{\mathbf{n}}{||\mathbf{n}||}$ where $\mathbf{n}$ is a white Gaussian vector. Clearly we have $\mathbf{u}^T\mathbf{u}=||\mathbf{u}||^2=1$. Moreover, for independent vectors $\mathbf{u}_1$ and $\mathbf{u}_2$, the distribution of $\mathbf{u}_1^T\mathbf{u}_2$ is known as discussed in here.

What is the distribution of $\mathbf{u}^T\mathbf{\Sigma}\mathbf{u}$ where $\mathbf{\Sigma}$ is a covariance matrix? This can also be seen as the norm of $\mathbf{\Sigma}^\frac{1}{2}\mathbf{u}$.

Thanks for any helps and hints.

Answer 1

I've been thinking about this question too and I haven't come across any definitive answer. One first glance, I doubt that the distribution of $Y=U^TAU=\frac{X^TAX}{X^TX},\; X\sim\mathcal{N}(0,I)$ has a simple closed form as even $X^TAX$ doesn't have one for general A, though it can be expressed by the characteristic function. However, for large enough dimension $n$, the distribution is approximately gaussian: let $X^TX = n\beta_X$, then $Y=\frac{1}{n}\sum_{i}\left(\sum_j\beta_x^{-1} A_{ij}x_ix_j\right)$ and the approximation follows from the central limit theorem). So from a practical standpoint, we might be able to estimate the distribution as $\mathcal{N}(E(Y),\text{var}(Y))$. Even so, I haven't come across an explicit expression for the mean/variance of $Y$. The best I can see is either 1) use lower/upper bounds or 2) approximate with Taylor's theorem. The first route yields for positive definite A: $$2\left(1-\frac{1}{\sqrt{n-4}}\right)\frac{\text{tr }A}{n-2}\leq E(Y) \leq 2\left(1+\frac{1}{\sqrt{n-4}}\right)\frac{\text{tr }A}{n-2}$$ using the fact that in general $E(X_1X_2) = \text{cov}(X_1,X_2) + E(X_1)E(X_2)$ and $|\text{cov}(X_1,X_2)| \leq \sqrt{\text{var}(X_1)\text{var}(X_2)}$ with $$\begin{equation}E(X^TAX) = \text{tr }A \\E([X^TX]^{-1})=\frac{1}{n-2}\\\text{var}(X^TAX)=2\text{ tr }A^2\\\text{var}([X^TX]^{-1})=\frac{2}{(n-2)^2(n-4)}\end{equation}$$ and the relation that $\text{tr }A^2 \leq (\text{tr }A)^2$ for SPD matrices. I haven't gotten around to it yet, but I believe similar bounds for the variance can be found using the expressions in this post, but will probably require approximating $X^TAX$ by some distribution that is easier to deal with unless A is a matrix of integers.

I actually just thought of another approach, to find the pdf of $Y^\prime = [Y,x_2,...,x_N]$ using the standard change of variable formula and then integrate over the unwanted $x_i$ to get the marginal distribution of $Y$. The resulting integral, however, doesn't look tractable.

EDIT: I recently returned to this problem for distributional bounds on the change to singular values of a perturbed matrix $A + \delta E$ where $E$ is diagonal and $\text{diag}(E)$ constrained to the nullspace of $(U_k \circ U_k)^T$. Anyway, since my original post I've had some dealings with elliptical distributions, of which this distribution is one example. Again exploiting the connection between normal and unit norm vectors, we have the $x = \rho u$ where $\rho ^2\sim\chi_k^2$. As @enthdegree links, the expected value is $\frac{1}{n}\text{tr }A$, and the variance is $$\frac{n \text{tr }A^2-(\text{tr }A)^2 }{n^2(1+\frac{1}{2}n)}$$ Higher (non-central) moments are also available by solving $E[\rho ^{2k}]E[(u^TAu)^k] = E[(x^TAx)^k]$, that is its non-central moments are related to those of the generalized chi-squared distribution by the non-central moments of the Chi-squared distribution. I found a paper that gives a moment/cumulant generating function for the product of quadratic forms that could be used (with some effort) for a moment matching approximation when the normal approximation isn't enough.

Answer 2

Let $P\Lambda P^\top$ be the eigendecomposition of $\mathbf\Sigma$. Then $$ R:=\mathbf{u}^T\mathbf{\Sigma}\mathbf{u}=\frac{Y^\top \Lambda Y}{Y^\top Y}=\frac{\sum_{i=1}^n\lambda_iY_i^2}{\sum_{i=1}^n Y_i^2}, $$ where $Y:=P^\top \mathbf{n}$, and $Y_i^2\overset{\text{i.i.d.}}{\sim}\chi_1^2$. It can be shown that $R$ is independent of $Y^{\top}Y$, which implies that $$ \mathsf{E}R^p=\frac{\mathsf{E}[Y^{\top}\Lambda Y]^p}{\mathsf{E}[Y^\top Y]^p}. $$ For example, $\mathsf{E}R=\frac{\operatorname{tr}(\mathbf{\Sigma})}{n}$, and $\operatorname{Var}(R)=\frac{2(n\operatorname{tr}(\mathbf{\Sigma}^2)-\operatorname{tr}(\mathbf{\Sigma})^2)}{n^2(n+2)}$.

---

As for the distribution of $R$, consider the eigendecomposition $\mathbf{\Sigma}-rI=P\Lambda_r P^\top$ and note that \begin{align} \mathsf{P}(R\le r)&=\mathsf{P}(\mathbf{n}^\top(\mathbf{\Sigma}-rI)\mathbf{n}\le 0)=\mathsf{P}(Y^\top \Lambda_r Y\le 0) \\ &=\mathsf{P}\!\left(\sum_{i=1}^n\lambda_{r,i} Y_i^2\le 0\right), \end{align} where $Y_i^2\overset{\text{i.i.d.}}{\sim}\chi_1^2$.