Draw n-dimensional uniform sample from a unit (n-1)-sphere defined by (n-1)-dimensional variable So my question is same as this and this which were not answered.
I will explain why do I need this. Perhaps, there could also be a different solution. I use a black box gradient-free optimizer (Bayesian Optimization for clarity) of some cost function $f(x)$, where $x \in R^n$. A special property of the function is that $f(x) = f(ax)$ for any scalar $ a>0$, so the relative scale of $x$ does not matter. Therefore, I add a constraint $||x||_2=1$ which well-poses the search problem.
Possible options are:
1) Black box optimizers usually avoid such constraints. But it is possible to add a regularizer like $(||x||_2-1)^2$ to the cost function. However, I think it is possible to avoid it using the third option.
2) Normalize all elements of $x$ as $x \gets x / ||x||_2$. The inverse mapping is one-to-many, which means that BO will waste time sampling $ax, a>0$ vectors, if for some $a$ we already know $f(ax)$.
For the same reason, this is not applicable in my case.
3) Reduce BO sampling dimension by 1, that is to sample $(n-1)$-dimensional variables, which on a unit (n-1)-sphere map into $n$-dimensional Cartesian vector. I am trying to figure out how to do it, but cannot find any way of doing it.
 A: Using spherical coordinates, 
\begin{align*}
x_1 &= \cos(\phi_1)\,,\\
x_2 &= \sin(\phi_1)\cos(\phi_2)\,,\\
&\quad\quad\vdots\\
x_{n-1} &= \sin(\phi_1)\cdots\sin(\phi_{n-2})\cos(\phi_{n-1})\,,\\
x_n&=\sin(\phi_1)\cdots\sin(\phi_{n-2})\sin(\phi_{n-1})\,,\\
\end{align*}
the $n$-dimensional unit sphere is parameterised by $(n-1)$ angles, $\phi_1,\ldots,\phi_{n-2}\in(0,\pi)$ and $\phi_{n-1}\in(0,2\pi)$. The uniform distribution on that sphere is associated with the joint density
$$\sin(\phi_1)^{n-2}\sin(\phi_2)^{n-3}\cdots\sin(\phi_{n-2})$$
which means that the $\phi_i$'s are independent with respective densities
$$\sin(\phi_1)^{n-2}\,,\ \sin(\phi_2)^{n-3}\,,\ldots,\ \sin(\phi_{n-2})$$and the Uniform $U(0,2\pi)$ on $\phi_{n-1}$. Simulating directly the $\phi_i$'s according to these densities thus returns the proper distribution of the angles for a uniform on the sphere, from which the Euclidean coordinates can be derived.

Simulating the $\phi_i$'s can be done (rather inefficiently) by accept-reject when using $\pi_1(\phi)\propto\sin(\phi)$ as a proposal since
$$\sin(\phi)^d \le \sin(\phi)\qquad d\ge 1\qquad \phi\in(0,\pi)$$as illustrated by the following picture (representing $\sin^d$ for $d=1,2,\cdots,5$ and simulations from $\pi_1$ accepted or rejected depending on the target: 

There exists however a more efficient way to simulate from the densities $\sin(\phi)^d$. Indeed, if$$\phi\sim\pi_d(\phi)\propto\sin(\phi)^d$$then, by a change of variable$$X=\cos(\phi)\sim f_d(x)\propto (1-x^2)^{(d-1)/2}\mathbb{I}_{(-1,1)}(x)$$and$$Y=X^2\sim g_d(y)\propto (1-y)^{(d-1)/2}y^{-1/2}\mathbb{I}_{(0,1)}(x)$$meaning$$Y\sim\mathcal{B}e(1/2,(d+1)/2)$$Therefore simulating $\phi$ from $\pi_d$ amounts to


*

*generate $Y\sim\mathcal{B}e(1/2,(d+1)/2)$ and $S=1-2\mathbb{I}_{U<1/2}$ for $U\sim\mathcal{U}(0,1)$

*take $\phi=\cos^{-1}(S*\sqrt{Y})$
As shown by the picture below for $d=3$, the fit on $10^5$ simulations is correct, as expected:

