Can anybody suggest how I can compute the second moment (or the whole moment generating function) of the cosine of two gaussian random vectors $x,y$, each distributed as $\mathcal N (0,\Sigma)$, independent of each other? IE, moment for the following random variable

$$\frac{\langle x, y\rangle}{\|x\|\|y\|}$$

The closest question is Moment generating function of the inner product of two gaussian random vectors which derives MGF for the inner product. There's also this answer from mathoverflow which links this question to distribution of eigenvalues of sample covariance matrices, but I don't immediately see how to use those to compute the second moment.

I suspect that second moment scales in proportion to half-norm of eigenvalues of $\Sigma$ since I get this result through algebraic manipulation for 2 dimensions, and also for 3 dimensions from guess-and-check. For eigenvalues $a,b,c$ adding up to 1, second moment is:


Using the following for numerical check

val1[a_, b_, c_] := (a + b + c)/(Sqrt[a] + Sqrt[b] + Sqrt[c])^2
val2[a_, b_, c_] := Block[{},
  x := {x1, x2, x3};
  y := {y1, y2, y3};
  normal := MultinormalDistribution[{0, 0, 0}, ( {
      {a, 0, 0},
      {0, b, 0},
      {0, 0, c}
     } )];
  vars := {x \[Distributed] normal, y \[Distributed] normal};
  NExpectation[(x.y/(Norm[x] Norm[y]))^2, vars]]

  val1[1.5,2.5,3.5] - val2[1.5,2.5,3.5]

Checking the formula for 4 variables (within numerical bounds):

val1[a_, b_, c_, 
  d_] := (a + b + c + d)/(Sqrt[a] + Sqrt[b] + Sqrt[c] + Sqrt[d])^2
val2[a_, b_, c_, d_] := Block[{},
  x := {x1, x2, x3, x4};
  y := {y1, y2, y3, y4};
  normal := 
   MultinormalDistribution[{0, 0, 0, 
     0}, {{a, 0, 0, 0}, {0, b, 0, 0}, {0, 0, c, 0}, {0, 0, 0, d}}];
  vars := {x \[Distributed] normal, y \[Distributed] normal};
  NExpectation[(x.y/(Norm[x] Norm[y]))^2, vars]]

val1[0.5, 1.5, 2.5, 3.5] - val2[0.5, 1.5, 2.5, 3.5]
  • $\begingroup$ Due to rotational freedom, since the cosine is invariant under rotations, one of the vectors can be assumed to be a unit vector in whatever direction is most convenient. That should simplify the problem quite a bit, to the second moment of the cosine of $x \in \mathcal{N}(0,\Sigma)$ with respect to $(1,0,0,\ldots)$. EDIT: Actually, this depends on the symmetry of $\Sigma$. $\endgroup$
    – jwimberley
    Commented Feb 28, 2017 at 22:32
  • 2
    $\begingroup$ w huber's answer here may be of interest: stats.stackexchange.com/a/85977/37483 $\endgroup$
    – KOE
    Commented Mar 1, 2017 at 2:57
  • $\begingroup$ @Student001 indeed, the 1/n rate derived in that question seems to be a special case of this formula, since we remove a degree of freedom by normalizing trace of covariance matrix to 1 $\endgroup$ Commented Mar 2, 2017 at 0:50
  • $\begingroup$ Aside: Note that, wlog, $\Sigma$ is diagonal. $\endgroup$
    – cardinal
    Commented Mar 3, 2017 at 6:46
  • $\begingroup$ I found the question of the distribution of $\frac{x}{\|x\|}$ being asked at least 3 times on crossvalidated, so hopefully this post will popularize the notion of "projected normal distribution" so it is no longer a question! :) $\endgroup$
    – Henry.L
    Commented Mar 7, 2017 at 3:25

1 Answer 1


Hey Yaroslav, you really do not have to hurry accepting my answer on MO and are more than welcomed to ask further details :).

Since you reformulate the question in 3-dim I can see exactly what you want to do. In MO post I thought you only need to calculate the largest cosine between two random variables. Now the problem seems tougher.

First, we calculate the normalized Gaussian $\frac{X}{\|X\|}$, which is not a trivial job since it actually has a name "projected normal distribution" because we can rewrite the multivariate normal density $X$ in terms of its polar coordinate $(\|X\|,\frac{X}{\|X\|})=(r,\boldsymbol{\theta})$. And the marginal density for $\boldsymbol{\theta}$ can be obtained in $$\int_{\mathbb{R}^{+}}f(r,\boldsymbol{\theta})dr$$

An important instance is that in which $x$ has a bivariate normal distribution $N_2(\mu,\Sigma)$, in which $\|x\|^{-1}x$ is said to have a projected normal (or angular Gaussian or offset normal) distribution.[Mardia&Peter]p.46

In this step we can obtain distributions $\mathcal{PN}_{k}$ for $\frac{X}{\|X\|}\perp\frac{Y}{\|Y\|}$, and hence their joint density $(\frac{X}{\|X\|},\frac{Y}{\|Y\|})$ due to independence. As for a concrete density function of projected normal distribution, see [Mardia&Peter] Chap 10. or [2] Equation (4) or [1] . (Notice that in [2] they also assume a special form of covariance matrix $\Sigma=\left(\begin{array}{cc} \Gamma & \gamma\\ \gamma' & 1 \end{array}\right)$)

Second, since we already obtained their joint density, their inner product can be readily derived using transformation formula $$\left(\frac{X}{\|X\|},\frac{Y}{\|Y\|}\right)\mapsto\frac{X}{\|X\|}\cdot\frac{Y}{\|Y\|}$$. Also see [3].

As long as we computed the density, the second moment is only a problem of integration.


[Mardia&Peter]Mardia, Kanti V., and Peter E. Jupp. Directional statistics. Vol. 494. John Wiley & Sons, 2009.

[1]Wang, Fangpo, and Alan E. Gelfand. "Directional data analysis under the general projected normal distribution." Statistical methodology 10.1 (2013): 113-127.

[2]Hernandez-Stumpfhauser, Daniel, F. Jay Breidt, and Mark J. van der Woerd. "The general projected normal distribution of arbitrary dimension: modeling and Bayesian inference." Bayesian Analysis (2016). https://projecteuclid.org/download/pdfview_1/euclid.ba/1453211962

[3]Moment generating function of the inner product of two gaussian random vectors

  • $\begingroup$ @YaroslavBulatov Hopefully this is well worth your bounty! $\endgroup$
    – Henry.L
    Commented Mar 7, 2017 at 3:27
  • $\begingroup$ The answer I posted on MO is not exactly what the OP wanted because I was thinking that he is searching for the canonical angle. my bad. $\endgroup$
    – Henry.L
    Commented Mar 7, 2017 at 3:29
  • 1
    $\begingroup$ Could you provide a proof that assuming identity covariance matrix is w.l.o.g? It's not obvious to me. It's "easy" to show cardinal's claim that diagonal matrix is w.l.o.g, but how do you get rid of the eigenvalues? $\endgroup$
    – KOE
    Commented Mar 7, 2017 at 14:13
  • $\begingroup$ @Student001 If $\Sigma=P'\Lambda P$, then $PX$ have an identity covariance matrix. $\endgroup$
    – Henry.L
    Commented Mar 7, 2017 at 16:26
  • 1
    $\begingroup$ No, if $P'\Lambda P$ is the spectral decomposition of $\Sigma$, then $PX$ as covariance matrix $\Lambda$, which need not be the identity, so at least that step doesn't justify $\Sigma = I$ w.l.o.g. Maybe your last comment does, I'm not sure. $\endgroup$
    – KOE
    Commented Mar 7, 2017 at 16:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.