Note on the normalization constant

In the textbook, Directional Statistics ^[3] by Mardia and Jupp, the normalization constant given for the Von Mises Fisher probability density is apparently different from the one given here: ${\displaystyle C_{p}(\kappa )}$. In that book, the normalization constant is specified as:

${\displaystyle C_{p}^{*}(\kappa )={\frac {({\frac {\kappa }{2}})^{p/2-1}}{\Gamma (p/2)I_{p/2-1}(\kappa )}}}$

where ${\displaystyle \Gamma }$ is the gamma function. This is resolved by noting that Mardia and Jupp give the density "with respect to the uniform distribution", while the density here is specified in the usual way, with respect to Lebesgue measure. The density (w.r.t. Lebesgue measure) of the uniform distribution is the reciprocal of the surface area of the (p-1)-sphere, so that the uniform density function is given by the constant:

${\displaystyle C_{p}(0)={\frac {\Gamma (p/2)}{2\pi ^{p/2}}}}$

It then follows that:

${\displaystyle C_{p}^{*}(\kappa )={\frac {C_{p}(\kappa )}{C_{p}(0)}}}$

While the value for ${\displaystyle C_{p}(0)}$ was derived above via the surface area, the same result may be obtained by setting ${\displaystyle \kappa =0}$ in the above formula for ${\displaystyle C_{p}(\kappa )}$. This can be done by noting that the series expansion for ${\displaystyle I_{p/2-1}(\kappa )}$ divided by ${\displaystyle \kappa ^{p/2-1}}$ has but one non-zero term at ${\displaystyle \kappa =0}$. (To evaluate that term, one needs to use the definition ${\displaystyle 0^{0}=1}$.)

Support

The support of the Von Mises–Fisher distribution is the hypersphere, or more specifically, the ${\displaystyle (p-1)}$-sphere, denoted as

${\displaystyle S^{p-1}=\left\{\mathbf {x} \in \mathbb {R} ^{p}:\left\|\mathbf {x} \right\|=1\right\}}$

This is a ${\displaystyle (p-1)}$-dimensional manifold embedded in ${\displaystyle p}$-dimensional Euclidean space, ${\displaystyle \mathbb {R} ^{p}}$.

Mean direction

A series of N independent unit vectors ${\displaystyle x_{i}}$ are drawn from a von Mises–Fisher distribution. The maximum likelihood estimates of the mean direction ${\displaystyle \mu }$ is simply the normalized arithmetic mean, a sufficient statistic:^[3]

${\displaystyle \mu ={\bar {x}}/{\bar {R}},{\text{where }}{\bar {x}}={\frac {1}{N}}\sum _{i}^{N}x_{i},{\text{and }}{\bar {R}}=\|{\bar {x}}\|,}$

Concentration parameter

Use the modified Bessel function of the first kind to define

${\displaystyle A_{p}(\kappa )={\frac {I_{p/2}(\kappa )}{I_{p/2-1}(\kappa )}}.}$

Then:

${\displaystyle \kappa =A_{p}^{-1}({\bar {R}}).}$

Thus ${\displaystyle \kappa }$ is the solution to

${\displaystyle A_{p}(\kappa )={\frac {\left\|\sum _{i}^{N}x_{i}\right\|}{N}}={\bar {R}}.}$

A simple approximation to ${\displaystyle \kappa }$ is (Sra, 2011)

${\displaystyle {\hat {\kappa }}={\frac {{\bar {R}}(p-{\bar {R}}^{2})}{1-{\bar {R}}^{2}}},}$

A more accurate inversion can be obtained by iterating the Newton method a few times

${\displaystyle {\hat {\kappa }}_{1}={\hat {\kappa }}-{\frac {A_{p}({\hat {\kappa }})-{\bar {R}}}{1-A_{p}({\hat {\kappa }})^{2}-{\frac {p-1}{\hat {\kappa }}}A_{p}({\hat {\kappa }})}},}$

${\displaystyle {\hat {\kappa }}_{2}={\hat {\kappa }}_{1}-{\frac {A_{p}({\hat {\kappa }}_{1})-{\bar {R}}}{1-A_{p}({\hat {\kappa }}_{1})^{2}-{\frac {p-1}{{\hat {\kappa }}_{1}}}A_{p}({\hat {\kappa }}_{1})}}.}$

Standard error

For N ≥ 25, the estimated spherical standard error of the sample mean direction can be computed as:^[4]

${\displaystyle {\hat {\sigma }}=\left({\frac {d}{N{\bar {R}}^{2}}}\right)^{1/2}}$

where

${\displaystyle d=1-{\frac {1}{N}}\sum _{i}^{N}\left(\mu ^{T}x_{i}\right)^{2}}$

It is then possible to approximate a ${\displaystyle 100(1-\alpha )\%}$ a spherical confidence interval (a confidence cone) about ${\displaystyle \mu }$ with semi-vertical angle:

${\displaystyle q=\arcsin \left(e_{\alpha }^{1/2}{\hat {\sigma }}\right),}$ where ${\displaystyle e_{\alpha }=-\ln(\alpha ).}$

For example, for a 95% confidence cone, ${\displaystyle \alpha =0.05,e_{\alpha }=-\ln(0.05)=2.996,}$ and thus ${\displaystyle q=\arcsin(1.731{\hat {\sigma }}).}$

General case

An algorithm for drawing pseudo-random samples from the Von Mises Fisher (VMF) distribution was given by Ulrich^[5] and later corrected by Wood.^[6] An implementation in R is given by Hornik and Grün;^[7] and a fast Python implementation is described by Pinzón and Jung.^[8]

To simulate from a VMF distribution on the ${\displaystyle (p-1)}$-dimensional unitsphere, ${\displaystyle S^{p-1}}$, with mean direction ${\displaystyle {\boldsymbol {\mu }}\in S^{p-1}}$, these algorithms use the following radial-tangential decomposition for a point ${\displaystyle \mathbf {x} \in S^{p-1}\subset \mathbb {R} ^{p}}$ :

${\displaystyle \mathbf {x} =t{\boldsymbol {\mu }}+{\sqrt {1-t^{2}}}\mathbf {v} }$

where ${\displaystyle \mathbf {v} \in \mathbb {R} ^{p}}$ lives in the tangential ${\displaystyle (p-2)}$-dimensional unit-subsphere that is centered at and perpendicular to ${\displaystyle {\boldsymbol {\mu }}}$; while ${\displaystyle t\in [-1,1]}$. To draw a sample ${\displaystyle \mathbf {x} }$ from a VMF with parameters ${\displaystyle {\boldsymbol {\mu }}}$ and ${\displaystyle \kappa }$, ${\displaystyle \mathbf {v} }$ must be drawn from the uniform distribution on the tangential subsphere; and the radial component, ${\displaystyle t}$, must be drawn independently from the distribution with density:

${\displaystyle f_{\text{radial}}(t;\kappa ,p)={\frac {(\kappa /2)^{\nu }}{\Gamma ({\frac {1}{2}})\Gamma (\nu +{\frac {1}{2}})I_{\nu }(\kappa )}}e^{t\kappa }(1-t^{2})^{\nu -{\frac {1}{2}}}}$

where ${\displaystyle \nu ={\frac {p}{2}}-1}$. The normalization constant for this density may be verified by using:

${\displaystyle I_{\nu }(\kappa )={\frac {(\kappa /2)^{\nu }}{\Gamma ({\frac {1}{2}})\Gamma (\nu +{\frac {1}{2}})}}\int _{-1}^{1}e^{t\kappa }(1-t^{2})^{\nu -{\frac {1}{2}}}\,dt}$

as given in Appendix 1 (A.3) in Directional Statistics.^[3] Drawing the ${\displaystyle t}$ samples from this density by using a rejection sampling algorithm is explained in the above references. To draw the uniform ${\displaystyle \mathbf {v} }$ samples perpendicular to ${\displaystyle {\boldsymbol {\mu }}}$, see the algorithm in,^[8] or otherwise a Householder transform can be used as explained in Algorithm 1 in.^[9]

3-D sphere

To generate a Von Mises–Fisher distributed pseudo-random spherical 3-D unit vector^[10][11] ${\textstyle \mathbf {X} _{s}}$ on the ${\textstyle S^{2}}$sphere for a given ${\textstyle \mu }$ and ${\textstyle \kappa }$, define

${\displaystyle \mathbf {X} _{s}=[\phi ,\theta ,r]}$

where ${\textstyle \phi }$ is the polar angle, ${\textstyle \theta }$ the equatorial angle, and ${\textstyle r=1}$ the distance to the center of the sphere

for ${\textstyle \mathbf {\mu } =[0,(.),1]}$ the pseudo-random vector is then given by

${\displaystyle \mathbf {X} _{s}=[\arccos W,V,1]}$

where ${\textstyle V}$ is sampled from the continuous uniform distribution ${\textstyle U(a,b)}$ with lower bound ${\textstyle a}$ and upper bound ${\textstyle b}$

${\displaystyle V\sim U(0,2\pi )}$

and

${\displaystyle W=1+{\frac {1}{\kappa }}(\ln \xi +\ln(1-{\frac {\xi -1}{\xi }}e^{-2\kappa }))}$

where ${\textstyle \xi }$ is sampled from the standard continuous uniform distribution ${\textstyle U(0,1)}$

${\displaystyle \xi \sim U(0,1)}$

here, ${\textstyle W}$should be set to ${\textstyle W=1}$ when ${\textstyle \mathbf {\xi } =0}$ and ${\textstyle \mathbf {X} _{s}}$ rotated to match any other desired ${\textstyle \mu }$.

Component marginal of uniform distribution

For ${\displaystyle 1\leq i\leq p}$, let ${\displaystyle x_{i}}$ be any component of ${\displaystyle \mathbf {x} \in S^{p-1}}$. The marginal distribution for ${\displaystyle x_{i}}$ has the density:^[12][13]

${\displaystyle f_{i}(x_{i};p)=f_{\text{radial}}(x_{i};\kappa =0,p)={\frac {(1-x_{i}^{2})^{{\frac {p-1}{2}}-1}}{B{\bigl (}{\frac {1}{2}},{\frac {p-1}{2}}{\bigr )}}}}$

where ${\displaystyle B(\alpha ,\beta )}$ is the beta function. This distribution may be better understood by highlighting its relation to the beta distribution:

${\displaystyle {\begin{aligned}x_{i}^{2}&\sim {\text{Beta}}{\bigl (}{\frac {1}{2}},{\frac {p-1}{2}}{\bigr )}&&{\text{and}}&{\frac {x_{i}+1}{2}}&\sim {\text{Beta}}{\bigl (}{\frac {p-1}{2}},{\frac {p-1}{2}}{\bigr )}\end{aligned}}}$

where the Legendre duplication formula is useful to understand the relationships between the normalization constants of the various densities above.

Note that the components of ${\displaystyle \mathbf {x} \in S^{p-1}}$ are not independent, so that the uniform density is not the product of the marginal densities; and ${\displaystyle \mathbf {x} }$ cannot be assembled by independent sampling of the components.

Distribution of dot-products

In machine learning, especially in image classification, to-be-classified inputs (e.g. images) are often compared using cosine similarity, which is the dot product between intermediate representations in the form of unitvectors (termed embeddings). The dimensionality is typically high, with ${\displaystyle p}$ at least several hundreds. The deep neural networks that extract embeddings for classification should learn to spread the classes as far apart as possible and ideally this should give classes that are uniformly distributed on ${\displaystyle S^{p-1}}$.^[14] For a better statistical understanding of across-class cosine similarity, the distribution of dot-products between unitvectors independently sampled from the uniform distribution may be helpful.

Let ${\displaystyle \mathbf {x} ,\mathbf {y} \in S^{p-1}}$ be unitvectors in ${\displaystyle \mathbb {R} ^{p}}$, independently sampled from the uniform distribution. Define:

${\displaystyle {\begin{aligned}t&=\mathbf {x} '\mathbf {y} \in [-1,1],&r&={\frac {t+1}{2}}\in [0,1],&s&={\text{logit}}(r)=\log {\frac {1+t}{1-t}}\in \mathbb {R} \end{aligned}}}$

where ${\displaystyle t}$ is the dot-product and ${\displaystyle r,s}$ are transformed versions of it. Then the distribution for ${\displaystyle t}$ is the same as the marginal component distribution given above;^[13] the distribution for ${\displaystyle r}$ is symmetric beta and the distribution for ${\displaystyle s}$ is symmetric logistic-beta:

${\displaystyle {\begin{aligned}r&\sim {\text{Beta}}{\bigl (}{\frac {p-1}{2}},{\frac {p-1}{2}}{\bigr )},&s&\sim B_{\sigma }{\bigl (}{\frac {p-1}{2}},{\frac {p-1}{2}}{\bigr )}\end{aligned}}}$

The means and variances are:

${\displaystyle {\begin{aligned}E[t]&=0,&E[r]&={\frac {1}{2}},&E[s]&=0,\end{aligned}}}$

and

${\displaystyle {\begin{aligned}{\text{var}}[t]&={\frac {1}{p}},&{\text{var}}[r]&={\frac {1}{4p}},&{\text{var}}[s]&=2\psi '{\bigl (}{\frac {p-1}{2}}{\bigr )}\approx {\frac {4}{p-1}}\end{aligned}}}$

where ${\displaystyle \psi '=\psi ^{(1)}}$ is the first polygamma function. The variances decrease, the distributions of all three variables become more Gaussian, and the final approximation gets better as the dimensionality, ${\displaystyle p}$, is increased.

Matrix Von Mises-Fisher

The matrix von Mises-Fisher distribution (also known as matrix Langevin distribution^[15][16]) has the density

${\displaystyle f_{n,p}(\mathbf {X}$ ;\mathbf {F} )\propto \exp(\operatorname {tr} (\mathbf {F} ^{\mathsf {T}}\mathbf {X} ))}

supported on the Stiefel manifold of ${\displaystyle n\times p}$ orthonormal p-frames ${\displaystyle \mathbf {X} }$, where ${\displaystyle \mathbf {F} }$ is an arbitrary ${\displaystyle n\times p}$ real matrix.^[17][18]

Saw distributions

Ulrich,^[5] in designing an algorithm for sampling from the VMF distribution, makes use of a family of distributions named after and explored by John G. Saw.^[19] A Saw distribution is a distribution on the ${\displaystyle (p-1)}$-sphere, ${\displaystyle S^{p-1}}$, with modal vector ${\displaystyle {\boldsymbol {\mu }}\in S^{p-1}}$ and concentration ${\displaystyle \kappa \geq 0}$, and of which the density function has the form:

${\displaystyle f_{\text{Saw}}(\mathbf {x}$ ;{\boldsymbol {\mu }},\kappa )={\frac {g(\kappa \mathbf {x} '{\boldsymbol {\mu }})}{K_{p}(\kappa )}}}

where ${\displaystyle g}$ is a non-negative, increasing function; and where ${\displaystyle K_{P}(\kappa )}$ is the normalization constant. The above-mentioned radial-tangential decomposition generalizes to the Saw family and the radial compoment, ${\displaystyle t=\mathbf {x} '{\boldsymbol {\mu }}}$ has the density:

${\displaystyle f_{\text{Saw-radial}}(t;\kappa )={\frac {2\pi ^{p/2}}{\Gamma (p/2)}}{\frac {g(\kappa t)(1-t^{2})^{(p-3)/2}}{B{\bigl (}{\frac {1}{2}},{\frac {p-1}{2}}{\bigr )}K_{p}(\kappa )}}.}$

where ${\displaystyle B}$ is the beta function. Also notice that the left-hand factor of the radial density is the surface area of ${\displaystyle S^{p-1}}$.

By setting ${\displaystyle g(\kappa \mathbf {x} '{\boldsymbol {\mu }})=e^{\kappa \mathbf {x} '{\boldsymbol {\mu }}}}$, one recovers the VMF distribution.