Sampling on the sphere from $f(x) \propto x^TAx$

TL;DR

This paper introduces a sampling method for generating random unit vectors on a sphere according to a density proportional to quadratic forms, with an accompanying R implementation.

Contribution

It presents a novel sampling technique for vectors on the sphere based on quadratic form densities, including a practical R implementation.

Findings

Effective sampling from the specified distribution on the sphere.

Implementation in R facilitates practical application.

Potential use in statistical and machine learning models.

Abstract

A method for drawing random samples of unit vectors $x$ in $R^p$ with density proportional to $x^TAx$ where $A$ is a symmetric, positive definite matrix. Includes an R function which implements the method.

Equations

$$A=R\Lambda R^T$$

$$f(x|\kappa ,A)$$

$$f(u|\kappa ,\Lambda )=\frac{p}{\text{tr}(\Lambda )}{u}^T\Lambda u=\frac{p}{{\sum }_{j=1}^p{\lambda }_j}\sum _{j=1}^p{\lambda }_j{u}_j^2$$

$$u=\begin{bmatrix}u_{1}\\ u_{2}\\ \vdots\\ u_{j}\\ \vdots\\ u_{p-2}\\ u_{p-1}\\ u_{p}\end{bmatrix}=\begin{bmatrix}t_{1}\\ (1-t_{1}^{2})^{1/2}t_{2}\\ \vdots\\ \bigl{[}\prod_{\ell=1}^{j-1}(1-t_{\ell}^{2})^{1/2}\bigr{]}t_{j}\\ \vdots\\ \bigl{[}\prod_{\ell=1}^{p-3}(1-t_{\ell}^{2})^{1/2}\bigr{]}t_{p-2}\\ \bigl{[}\prod_{\ell=1}^{p-2}(1-t_{\ell}^{2})^{1/2}\bigr{]}\cos\phi_{p-1}\\ \bigl{[}\prod_{\ell=1}^{p-2}(1-t_{\ell}^{2})^{1/2}\bigr{]}\sin\phi_{p-1}\end{bmatrix}=\begin{bmatrix}t_{1}\\ b_{2}^{1/2}t_{2}\\ \vdots\\ b_{j}^{1/2}t_{j}\\ \vdots\\ b_{p-2}^{1/2}t_{p-2}\\ b_{p-1}^{1/2}\cos\phi_{p-1}\\ b_{p-1}^{1/2}\sin\phi_{p-1}\end{bmatrix}$$

$$b_1$$

$$b_j(t_1,\dots ,t_{j-1})$$

$$a_1$$

$$a_j(t_1,\dots ,t_{j-1})$$

$$f_{\text{Uniform}}(t_1,\dots ,t_{p-2},{\varphi }_{p-1})$$

$$g_k(t_1,\dots ,t_k)$$

$$f(t_1,\dots ,t_{p-2},{\varphi }_{p-1})$$

$$=$$

$$f(t_{1},\ldots,t_{p-2})=g_{p-2}(t_{1},\ldots,t_{p-2})\frac{p}{\sum_{k=1}^{p}\lambda_{k}}\left(a_{p-1}+\mbox{$\frac{1}{2}$}b_{p-1}(\lambda_{p-1}+\lambda_{p})\right)$$

$$f(t_1,\dots ,t_j)=g_j(t_1,\dots ,t_j)\frac{p}{{\sum }_{k=1}^p{\lambda }_k}\left(a_{j+1}+\frac{b_{j+1}}{p-j}\left(\sum _{k=j+1}^p{\lambda }_k\right)\right)$$

$$f(t_1)$$

$${\pi }_1=\frac{{\lambda }_1}{{\sum }_{j=1}^p{\lambda }_j}\text{and}{\pi }_2=\frac{{\sum }_{k=2}^p{\lambda }_k}{{\sum }_{k=1}^p{\lambda }_k}$$

$$B^{\ast}_{1}\sim\text{Beta}(\mbox{$\frac{3}{2}$},\mbox{$\frac{p-1}{2}$})$$

$$B^{\ast}_{1}\sim\text{Beta}(\mbox{$\frac{1}{2}$},\mbox{$\frac{p+1}{2}$})$$

$$f(t_j|t_1,\dots ,t_{j-1})$$

$${\pi }_{j1}$$

$${\pi }_{j2}$$

$${\pi }_{j3}$$

$$B^{\ast}_{j}\sim\text{Beta}(\mbox{$\frac{1}{2}$},\mbox{$\frac{p-j}{2}$})$$

$$B^{\ast}_{j}\sim\text{Beta}(\mbox{$\frac{3}{2}$},\mbox{$\frac{p-j}{2}$})$$

$$B^{\ast}_{j}\sim\text{Beta}(\mbox{$\frac{1}{2}$},\mbox{$\frac{p-j+2}{2}$})$$

$$f({\varphi }_{p-1}|t_1,\dots ,t_{p-2})$$

$$\displaystyle\mbox{$\frac{1}{4}$}U_{p-1}$$

$$=$$

Full text

Sampling on the sphere from $f(x)\propto x^{T}Ax$

Richard Arnold

School of Mathematics and Statistics

Victoria University of Wellington

Wellington, New Zealand

richard.arnold@vuw.ac.nz

1 Preliminaries

We wish to sample on the surface of the sphere $S^{p-1}=\{x:x\in\mathbb{R}^{p},\ x^{T}x=1\}$ from an axial distribution with density proportional to $x^{T}Ax$.

We assume that the matrix $A$ is symmetric, positive definite and can be diagonalised

[TABLE]

for orthogonal matrix $R$ and diagonal $\Lambda=\text{diag}(\lambda_{1},\ldots,\lambda_{p})$ with $\lambda_{j}>0$ for all $j=1,\ldots,p$.

Using the notation from (Arnold & Jupp, 2013) we can construct the density from an axial cardoid distribution

[TABLE]

where the last line follows from the requirement that $1-\kappa\text{tr}(A)/p=0$ which ensures that the desired proportionality holds.

In order to draw from this distribution we first take a draw $u$ from

[TABLE]

and then compute $x=Ru$.

2 Construction

The $p-$dimensional vector $u$ has $p-1$ degrees of freedom which we parameterise as $(t_{1},\ldots,t_{p-2},\phi_{p-1})$ with $t_{j}\in[-1,1]$ for $j=1,\ldots,p-2$ and $\phi_{p-1}\in[0,2\pi)$ (see, e.g. Mardia & Jupp, 2000, Ch. 9).

We then construct the vector $u$ as follows:

[TABLE]

Here we have defined the symbols

[TABLE]

and we further define

[TABLE]

Under this parameterisation the uniform distribution on the sphere has the form

[TABLE]

where $B(a,b)$ is the Beta function $\Gamma(a)\Gamma(b)/\Gamma(a+b)$ and we define

[TABLE]

If the density is given by (2) then with respect to our chosen parameterisation of $u$ we have

[TABLE]

Integrating with respect to $\phi_{p-1}\in[0,2\pi)$ we obtain the marginal distribution over $(t_{1},\ldots,t_{p-2})$:

[TABLE]

Further integrations over $t_{j}\in[-1,+1]$ lead to the general expression

[TABLE]

which holds for $j=1,\ldots,p-2$. When $j=1$ we ultimately have

[TABLE]

This latter distribution is a mixture of two Beta densities in $t_{1}^{2}$. To simulate from it we proceed as follows:

Compute the mixture probabilities

[TABLE] 2. 2.

Draw $X_{1}\sim\text{Categorical}((1,2);(\pi_{1},\pi_{2}))$ 3. 3.

If $X_{1}=1$ draw

[TABLE]

otherwise if $X_{1}=2$ draw

[TABLE] 4. 4.

Draw D_{1}\sim\text{Categorical}((-1,+1);(\mbox{\frac{1}{2}},\mbox{\frac{1}{2}})) 5. 5.

Set $t_{1}=D_{1}\sqrt{B^{\ast}_{1}}$.

To simulate $t_{2},\ldots,t_{p-2}$ we draw each $t_{j}$ conditional on $t_{1},\ldots,t_{j-1}$ by constructing the conditional distributions

[TABLE]

This is a mixture of three Beta densities in $t_{j}^{2}$. The mixing proportions are

[TABLE]

corresponding to the three component densities \text{Beta}(\mbox{\frac{1}{2}},\mbox{\frac{p-j}{2}}), \text{Beta}(\mbox{\frac{3}{2}},\mbox{\frac{p-j}{2}}) and \text{Beta}(\mbox{\frac{1}{2}},\mbox{\frac{p-j+2}{2}}) respectively.

Thus to simulate $t_{j}|t_{1},\ldots,t_{j-1}$

Compute the mixture probabilities $\pi_{j1}$, $\pi_{j2}$ and $\pi_{j3}$ using (17) 2. 2.

Draw $X_{j}\sim\text{Categorical}((1,2,3);(\pi_{j1},\pi_{j2},\pi_{j3}))$ 3. 3.

If $X_{j}=1$ draw

[TABLE]

otherwise if $X_{j}=2$ draw

[TABLE]

otherwise if $X_{j}=3$ draw

[TABLE] 4. 4.

Draw D_{j}\sim\text{Categorical}((-1,+1);(\mbox{\frac{1}{2}},\mbox{\frac{1}{2}})) 5. 5.

Set $t_{j}=D_{j}\sqrt{B^{\ast}_{j}}$.

We carry out this procedure for $j=2,\ldots,p-2$. Note that the procedure is also defined for $j=1$, in which case $\pi_{1}=0$ and we have the same two component mixture derived earlier.

The final step is to draw the angle $\phi_{p-1}\in[0,2\pi)$. Its conditional distribution is given by

[TABLE]

This density has four-fold symmetry: $\pm\phi$ and$\pi\pm\phi$ are all equivalent. To draw a value for $\phi_{p-1}$ we first draw $U_{p-1}\sim\text{Uniform}(0,1)$ and then solve the equation

[TABLE]

in the interval \phi^{\ast}\in[0,\mbox{\frac{\pi}{2}}]. We finally set $\phi_{p-1}$ to one of the four values $\{\pm\phi^{\ast},\pi\pm\phi^{\ast}\}$ with equal probability.

Once we have a full draw of the values $(t_{1},\ldots,t_{p-2},\phi_{p-1})$ we use (3) to construct the vector $u$, and finally set $x=Ru$ as the random draw from $f(x)\propto x^{T}Ax$.

3 R code

The following R function can be used to take $n={\tt n}$ draws from the density (1) where $A={\tt amat}$ is a $p\times p$ square, positive defnite matrix. The function returns an $n\times p$ matrix, the rows of which are the $n$ random vectors.

rxax <- function(n,amat) { # Eigen decomposition ee <- eigen(amat) dvec <- ee$values rmat <- ee$vectors csum <- rev(cumsum(rev(dvec))) p <- nrow(amat)

# Generate vectors in rotated frame where amat is diagonal
umat <- sapply(1:n,
               function(i) {
                  avec <- rep(NA,p-1)
                  bvec <- rep(NA,p-1)
                  tvec <- rep(NA,p-2)
                  # generate t[1]...t[p-2]
                  avec[1] <- 0
                  bvec[1] <- 1
                  for(j in 1:(p-2)) {
                      pivec <- c( avec[j]*(p-j+1),
                                  bvec[j]*dvec[j],
                                  bvec[j]*csum[j+1] )
                      pivec <- pivec/sum(pivec)
                      i <- sample(1:3, size=1, prob=pivec)
                      b <- rbeta(1,
                                 (c(1,3,1)/2)[i],
                                 (c(p-j,p-j,p-j+2)/2)[i])
                      d <- sample(c(-1,1), size=1)
                      tvec[j] <- d*b
                      bvec[j+1] <- bvec[j]*(1-tvec[j]^2)
                      avec[j+1] <- avec[j] + dvec[j]*bvec[j]*tvec[j]^2
                  }
                  # generate phi[p-1]
                  ustar <- runif(1)
                  ffunc <- function(phi, uval, a, b, d1, d2) {
                      c1 <- (b/4)*(d1-d2)
                      c2 <- a + 0.5*b*(d1+d2)
                      (1/(2*pi)*( phi + c1*sin(2*phi)/c2 ) - uval/4)
                  }
                  # generate angle in (0,pi/2)
                  phistar <- uniroot(ffunc, c(0,pi/2),
                                     uval=ustar, a=avec[p-1], b=bvec[p-1],
                                     d1=dvec[p-1], d2=dvec[p])$root
                  # randomise over fourfold symmetry
                  phistar <- (sample(c(0,1),size=1)*pi
                             +sample(c(-1,+1),size=1)*phistar)

                  # construct uvec
                  uvec <- c( sqrt(bvec)[1:(p-2)]*tvec,
                             sqrt(bvec[p-1])*c(cos(phistar), sin(phistar)) )
                  })
# rotate back to natural frame
xmat <- t(rmat%*%umat)
return(xmat)

}

An example use of the function is as follows: we construct a suitable symmetric $3\times 3$ matrix $A$ at random and then take 10 draws from it.

amat <- array(rnorm(9),dim=c(3,3)) amat <- t(amat)%*%amat xmat <- rxax(10, amat) xmat

         [,1]        [,2]       [,3]

[1,] 0.98674959 0.01738693 -0.1613163 [2,] -0.38065170 0.23323003 0.8948229 [3,] 0.24385388 -0.60136593 0.7608510 [4,] -0.44503351 -0.66007152 0.6051865 [5,] -0.85847885 -0.04198661 -0.5111274 [6,] 0.33440170 0.04373660 0.9414152 [7,] -0.09637544 -0.11457369 -0.9887288 [8,] 0.27959151 -0.70197708 0.6550242 [9,] 0.61130650 -0.18851719 0.7686128 [10,] -0.38150181 -0.76106892 0.5246241