Hyperspherical parameterization of the simplex

[sethaxen]

As noted in an earlier e-mail, if $z$ is a unit vector of length $N$, then squaring each element so $x_i = z_i^2$ gives us a point $x$ on the simplex. If all elements of $z$ are positive (i.e. $z$ is on the positive orthant of the hypersphere), then this map is bijective. So if we can generate unit vectors with positive elements, we can use that transform to also sample on the simplex.

As mentioned in #1, using exp to make each element positive in the end gives us the augmented softmax parameterization.

We can do this as well with the hyperspherical transform, which is written

$$ z_i = \left(\prod_{j=1}^{i-1} \sin(\phi_j)\right) \begin{cases} 1 & \text{ if } i = N\\ \cos(\phi_i) & \text{otherwise} \end{cases}, $$

where $\phi_j$ is in $(0, 2\pi]$ when $j=N-1$ and $(0, \pi]$ otherwise.

The requirement that $z_i$ be positive for $i < N-1$ is satisfied by requiring $\phi_i \in (0, \pi/2]$ for all $i$.

It turns out that the Jacobian determinant here is

$$ \begin{aligned} |J| &= 2^{N-1} \left(\prod_{i=1}^{N-2} \sin(\phi_i)^{N-i-1}\right) \left(\prod_{i=1}^N z_i \right)\\ &= 2^{N-1} \prod_{i=1}^{N-1}\tan(\phi_i) x_i\\ &= 2^{N-1} \prod_{i=1}^{N-1}\sin^{2(N-i)-1}(\phi_i) \cos(\phi_i). \end{aligned} $$

I suspect due to the trigonometric functions, it may not be extremely efficient. However, all of these trigonometric functions are bijective over this range, so we can replace sines or cosines with numbers constrained to the interval $(0, 1)$ to get a new transform.

I'll shortly add Stan functions for both of these.

[sethaxen]

Continuing further, if we let $s_i = \sin^2(\phi_i) \in (0, 1]$, then the transform from $s$ is

$$ x_i = \left(\prod_{j=1}^{i-1} s_j \right) \begin{cases} 1 & \text{ if } i = N\\ 1 - s_i & \text{otherwise} \end{cases}, $$

and its absolute Jacobian determinant is

$$ |J| = \prod_{i=1}^{N-2}s_i^{N-i - 1}. $$

Finally, letting $s_i = \operatorname{logistic}(y_i)$, we have the transform

$$ x_i = \left(\prod_{j=1}^{i-1} \operatorname{logistic}(y_j) \right) \begin{cases} 1 & \text{ if } i = N\\ 1 - \operatorname{logistic}(y_i) & \text{otherwise} \end{cases}, $$

$$ |J| = \prod_{i=1}^{N-1}\operatorname{logistic}(y_i)^{N-i} (1-\operatorname{logistic}(y_i)) $$

[mjhajharia]

we could reference https://arxiv.org/pdf/1010.3436.pdf

[sethaxen]

we could reference https://arxiv.org/pdf/1010.3436.pdf

Ah, perfect! Yes, the transformations are identical.

[mjhajharia]

yes!

sidenote: also because of some squashing that I did and the time when i wasn't merging PRs there were conflicts on all branches so I merged all the PRs and am fixing up that stuff now, if after an hour or so something's wrong let me know so I can figure out if I messed up something

[sethaxen]

Okay, thanks for the heads up!