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



 
 
Sampling on the sphere from f(x)∝ x^TAx
    Richard Arnold
 
        School of Mathematics and Statistics

        Victoria University of Wellington

        Wellington, New Zealand
 
        richard.arnold@vuw.ac.nz
    December 30, 2023
================================================================================================================================================================================




§ PRELIMINARIES


We wish to sample on the surface of the sphere S^p-1 = {x:
x∈ℝ^p,  x^Tx=1} from an axial distribution with density
proportional to x^TAx.

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

    A=RΛ R^T

for orthogonal matrix R and diagonal
Λ=diag(λ_1,…,λ_p) with λ_j>0 for all
j=1,…,p.

Using the notation from <cit.> we can construct the density
from an axial cardoid distribution

    f(x|κ,A)    =    1 + κ⟨ A, t([U])⟩
       =    1 + κtr(A(xx^T-I/p))
       =    1 + κ(x^TAx-tr(A)/p)
       =   p/tr(A) x^TAx

where the last line follows from the requirement that
1-κ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

    f(u|κ,Λ) = p/tr(Λ) u^TΛ u
                           = p/∑_j=1^p λ_j∑_j=1^p λ_ju_j^2

and then compute x=Ru.  




§ CONSTRUCTION


The p-dimensional vector u has p-1 degrees of freedom which we
parameterise as (t_1,…,t_p-2,ϕ_p-1) with t_j∈[-1,1] for
j=1,…,p-2 and ϕ_p-1∈[0,2π) <cit.>.

We then construct the vector u as follows:

    u = [   u_1;   u_2;     ⋮;   u_j;     ⋮; u_p-2; u_p-1;   u_p ]
         = [                               t_1;                 (1-t_1^2)^1/2 t_2;                                 ⋮;      [∏_ℓ=1^j-1(1-t_ℓ^2)^1/2] t_j;                                 ⋮;    [∏_ℓ=1^p-3(1-t_ℓ^2)^1/2] t_p-2; [∏_ℓ=1^p-2(1-t_ℓ^2)^1/2] cosϕ_p-1; [∏_ℓ=1^p-2(1-t_ℓ^2)^1/2] sinϕ_p-1 ]
         = [               t_1;       b_2^1/2 t_2;                 ⋮;       b_j^1/2 t_j;                 ⋮;   b_p-2^1/2 t_p-2; b_p-1^1/2cosϕ_p-1; b_p-1^1/2sinϕ_p-1 ]

Here we have defined the symbols

    b_1    =    1 
    
      b_j(t_1,…,t_j-1)    =   ∏_ℓ=1^j-1(1-t_ℓ^2)
             for j=2,…,p-1
       =    b_j-1(t_1,…,t_j-2) (1-t_j-1)^2

and we further define

    a_1    =    0 
    
      a_j(t_1,…,t_j-1)    =   ∑_k=1^j-1λ_k t_k^2
                                 ∏_ℓ=1^k-1(1-t_ℓ^2) 
       =   ∑_k=1^j-1λ_k b_k t_k^2
             for j=2,…,p-1
       =    a_j-1(t_1,…,t_j-2) + λ_j-1b_j-1t_j-1^2

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

    f_Uniform(t_1,…,t_p-2,ϕ_p-1)
         =   1/2π∏_j=1^p-2 B(12, p-j2)^-1
                     (1-t_j^2)^(p-j-2)/2
       =   1/2π g_p-2(t_1,…, t_p-2)

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

    g_k(t_1,…, t_k)    =   ∏_j=1^k B(12, p-j2)^-1
                     (1-t_j^2)^(p-j-2)/2    for k=1,…,p-2
       =    g_k-1(t_1,…, t_k-1) 
                     B(12, p-k2)^-1
                     (1-t_k^2)^(p-j-2)/2

If the density is given by (<ref>) then with respect to our
chosen parameterisation of u we have

    f(t_1,…,t_p-2,ϕ_p-1)      
       =   1/2π g_p-2(t_1,…, t_p-2) p/∑_j=1^pλ_j∑_j=1^p λ_j x_j^2
       =   1/2π g_p-2(t_1,…, t_p-2) p/∑_j=1^pλ_j(∑_j=1^p-2λ_j b_jt_j^2
                + 
                λ_p-1 b_p-1cos^2ϕ_p-1
                + 
                λ_p b_p-1sin^2ϕ_p-1)
       =   1/2π g_p-2(t_1,…, t_p-2) p/∑_j=1^pλ_j(a_p-1
                + 
                λ_p-1 b_p-1cos^2ϕ_p-1
                + 
                λ_p b_p-1sin^2ϕ_p-1)

Integrating with respect to ϕ_p-1∈[0,2π) we obtain the marginal
  distribution over (t_1,…,t_p-2):

    f(t_1,…,t_p-2) = 
      g_p-2(t_1,…, t_p-2) p/∑_k=1^pλ_k(a_p-1
                + 
                12 b_p-1(λ_p-1+λ_p)
          )

Further integrations over t_j∈[-1,+1] lead to the general expression

    f(t_1,…,t_j) = 
      g_j(t_1,…, t_j) p/∑_k=1^pλ_k(a_j+1
                + 
                b_j+1/p-j(∑_k=j+1^p λ_k)
          )

which holds for j=1,…,p-2.  When j=1 we ultimately have 

    f(t_1) 
         =    g_1(t_1) p/∑_k=1^pλ_k(a_2
                + 
                b_2/p-1(∑_k=2^p λ_k)
          )
       =    g_1(t_1) p/∑_j=1^pλ_j(λ_1t_1^2 + 
                1p-1(∑_k=2^p λ_k) (1-t_1^2))
       =   λ_1/∑_j=1^pλ_j
          B(32,p-12) t_1^2 (1-t_1^2)^(p-3)/2
         +∑_k=2^pλ_k/∑_k=1^pλ_k
          B(12,p+12) (1-t_1^2)^(p-1)/2

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
   
    π_1 = λ_1/∑_j=1^pλ_j    and    π_2 = ∑_k=2^pλ_k/∑_k=1^pλ_k


  * Draw X_1∼Categorical((1,2); (π_1,π_2))

  * If X_1=1 draw 
   
    B^∗_1 ∼Beta(32,p-12)

   otherwise if X_1=2 draw
   
    B^∗_1 ∼Beta(12,p+12)


  * Draw D_1∼Categorical((-1,+1); (12,12))

  * Set t_1=D_1√(B^∗_1).

To simulate t_2,…,t_p-2 we draw each t_j conditional on
t_1,…,t_j-1 by constructing the conditional distributions

    f(t_j|t_1,…,t_j-1)      
       =   f(t_1,…,t_j-1,t_j)/f(t_1,…,t_j-1)
       =   g_j(t_1,…, t_j)/g_j-1(t_1,…, t_j-1)
                a_j+1
                + 
                b_j+1/p-j(∑_k=j+1^p λ_k)
            /
                a_j
                + 
                b_j/p-j+1(∑_k=j^p λ_k)
            
       =   [(1-t_j^2)^(p+j-2)/2/B(12,p-j2)]
           [
           a_j + λ_jb_jt_j^2 
                     + 1p-j(∑_k=j+1^pλ_k)
                       b_j(1-t_j^2)
                 /
                 a_j + 1p-j+1(∑_k=j^pλ_k)
                       b_j
                ]

This is a mixture of three Beta densities in t_j^2.  The mixing proportions
are

    π_j1   =   a_j(p-j+1)/a_j(p-s+1) + b_j(∑_k=j^pλ_k)
    π_j2   =   b_jλ_j/a_j(p-s+1) + b_j(∑_k=j^pλ_k)
    π_j3   =   b_j(∑_k=j+1^pλ_k)/a_j(p-s+1) + b_j(∑_k=j^pλ_k)

corresponding to the three component densities 
Beta(12,p-j2), 
Beta(32,p-j2) and
Beta(12,p-j+22) respectively. 

Thus to simulate t_j|t_1,…,t_j-1


  * Compute the mixture probabilities π_j1, π_j2 and π_j3 using
   (<ref>)

  * Draw X_j∼Categorical((1,2,3); (π_j1,π_j2,π_j3))

  * If X_j=1 draw 
   
    B^∗_j ∼Beta(12,p-j2)

   otherwise if X_j=2 draw 
   
    B^∗_j ∼Beta(32,p-j2)

   otherwise if X_j=3 draw
   
    B^∗_j ∼Beta(12,p-j+22)


  * Draw D_j∼Categorical((-1,+1); (12,12))

  * Set t_j=D_j√(B^∗_j).

We carry out this procedure for j=2,…,p-2.  Note that the procedure is
also defined for j=1, in which case π_1=0 and we have the same two
component mixture derived earlier.

The final step is to draw the angle ϕ_p-1∈[0,2π).  Its conditional
  distribution is given by 

    f(ϕ_p-1|t_1,…,t_p-2)      
       =   f(t_1,…,t_p-2,ϕ_p-1)/f(t_1,…,t_p-2)
       =   1/2π
             a_p-1
                + 
                λ_p-1 b_p-1cos^2ϕ_p-1
                + 
                λ_p b_p-1sin^2ϕ_p-1/
             a_p-1
                + 
                12 b_p-1(λ_p-1+λ_p)
             
       =   1/2π[1 + 12b_p-1(λ_p-1 -λ_p)cos 2ϕ_p-1/
             a_p-1
                + 
                12 b_p-1(λ_p-1+λ_p)
             ]

This density has four-fold symmetry: ±ϕ andπ±ϕ are all
equivalent.  To draw a value for ϕ_p-1 we first draw
U_p-1∼Uniform(0,1) and then solve the equation

    14 U_p-1   =   ∫_0^ϕ^∗ f(ϕ'|t_1,…,t_p-2)  dϕ'
       =   1/2π[ϕ^∗ + 14b_p-1(λ_p-1 -λ_p)sin 2ϕ^∗/
             a_p-1
                + 
                12 b_p-1(λ_p-1+λ_p)
             ]
 
in the interval ϕ^∗∈[0,π2].  We finally set ϕ_p-1
to one of the four values {±ϕ^∗, π±ϕ^∗} with equal
probability.

Once we have a full draw of the values (t_1,…,t_p-2,ϕ_p-1) we 
use (<ref>) to construct the vector u, and finally set x=Ru
as the random draw from f(x)∝ x^TAx.  




§ R CODE


The following R function can be used to take n= n draws from the density
(<ref>) where A= amat is a p× p square, positive
defnite matrix.  The function returns an n× p matrix, the rows of which
are the n random vectors. 







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









agsm


