I'm working on an algorithm that is minimizing some loss function restricted to a polytope in $\mathbb{R}^n$:
$$ \mathcal{X}= \{ x \in \mathbb{R}^n : \sum_{i=1}^n x_i \leq 1, x_i \geq 0, \ i=1,\dots,n\}. $$
For a gradient estimating procedure, I want to sample from the cone of feasible directions $\Delta$ given any point $x \in \mathcal{X}$:
$$ \Delta = \{ \delta \in \mathbb{R}^n : x + \epsilon\delta \in \mathcal{X}, \|\delta\|_2 = 1, \text{for some} \ \epsilon > 0 \}. $$
This is trivial for all $x \in \text{int}(\mathcal{X})$, as we can simply sample from the $(n-1)$ unit sphere using standard normals $X_i, \text{for} \ i=1,\dots,n,$ and letting $\delta_i = X_i / \|X\|_2$. However, this becomes infeasible once $x$ is on the boundary of $\mathcal{X}$ and thus we have active constraints. Note that for active constraints $x_i = 0$ we can let $\delta_i = |X_i| / \|X\|_2$. However, when $\sum_{i=1}^n x_i = 1$, i.e., the constraint is active, this not sufficient. Applying an acceptance-rejection algorithm using the suggestion above will become extremely inefficient for $n$ large.
My question:
How to uniformly sample from $\Delta$ in case $\sum_{i=1}^n x_i = 1$?
The only helpful thing I could find was a paper on sampling from spherical triangles: https://www.graphics.cornell.edu/pubs/1995/Arv95c.pdf. I obviously need an algorithm for sampling from higher dimensional spherical polytopes.