While thinking of some algorithms related to machine learning, I went on a tangent and eventually asked myself if I could transform a standard normal distribution into a uniform distribution on the unit $n$-sphere. I originally thought this would be relatively easy, but it turned out to be more challenging than I originally thought.
The specific challenge I set for myself is as follows:
Suppose $X \in \mathbb{R}^n$, $X$ ~ $\mathcal{N}(\mathbf{0}, I)\;$ (meaning $X$ is an n-dimensional standard normal vector)
Find $f \;.\; \mathbb{R}^n \rightarrow \mathbb{R}^{n + 1}$ such that:
$f$ is continuous. Its domain is the $n$-dimensional Euclidean plane and its range is the unit $n$-sphere minus the point $-\mathbf{e}_{n + 1}$.
$f(X)$ is uniformly distributed on the $n$-sphere (minus $-\mathbf{e}_{n + 1}$)
$f$ is isotropic. This means for all $d > 0$, the set $S = \{f(x) \; | \; ||x|| = d\}$ remains constant in the $(n + 1)$th dimension.
$f(\mathbf{0}) = \mathbf{e}_{n + 1}$
When I search for similar topics online, I only find pages about the fact that an $(n + 1)$-dimensional standard normal variable divided by its norm yields a variable uniformly distributed on the $n$-sphere. This question is different because the input variable is itself $n$-dimensional, and for my application I'm not allowed to alter the data by appending an artificial standard normally distributed variable to forcibly make $(n + 1)$-dimensions.
Since the question relates to a machine learning algorithm I'm trying to create, ideally $f$ should have a closed form expression that is fast to numerically compute even with rather large $n$.
I tried to solve this for a long time and I'm stuck now. Does anyone have good ideas for how this function can be constructed? If not, is there a good way to quickly numerically estimate a function that can do the job?
What I tried so far:
The approach I thought had the most promise was to do an inverse stereographic projection. This is inspired by the fact that the 1D stereographic projection maps a uniform distribution on the unit circle to a Cauchy distribution on $\mathbb{R}$. This means I can transform normally distributed data into Cauchy distributed data and run the inverse stereographic projection.
However, I can't seem to find a simple way to describe the distribution for a stereographic projections of uniform distributions on higher dimensional spheres. It's not a multivariate Cauchy distribution. After some algebra, I seem to end up with some messy formulas that I can't make sense of. This is how I approached the problem:
Define polar coordinates for the $n$-sphere: $\theta_1 \in (-\pi, \pi], \; \theta_2, ...\theta_{n - 1}, \phi \in [0, \pi]$ where $\mathbf{y} \in \mathbb{R}^{n + 1} = (\sin(\theta_1)\sin(\theta_2)...\sin(\theta_{n - 1})\sin(\phi), \; \cos(\theta_1)\sin(\theta_2)...\sin(\theta_{n - 1})\sin(\phi), \; \cos(\theta_2)\sin(\theta_3)...\sin(\theta_{n - 1})\sin(\phi), \; ... \cos(\phi))$
Define polar coordinates for the $n$-dimensional plane: $\theta_1 \in (-\pi, \pi], \; \theta_2, ...\theta_{n - 1} \in [0, \pi], r \in \mathbb{R+}$ where $\mathbf{x} \in \mathbb{R}^n = (r \sin(\theta_1)\sin(\theta_2)...\sin(\theta_{n - 1}), \; ... r \cos(\theta_{n - 1}))$
$dA_{\text{sphere}} = \sin(\phi)^{n - 1} dL \, d\phi = $ The differential surface area element on the $n$-sphere, where $dL$ represents the differential surface area element of the "great" $(n - 1)$-sphere corresponding to $\phi = \pi / 2$
$dA_{\text{plane}} = r^{d - 1} dL \, dr = $ The differential surface area element on the plane.
Stereographic projection of the $n$-sphere from the point $-\mathbf{e}_{n + 1}$ to $\{(v, 0) \; | \; v \in \mathbb{R}^n\}$ preserves $\theta_1, ... \theta_{n - 1}$ and sets $r = \tan(\phi / 2)$, meaning $\phi = 2 \tan^{-1}(r)$
$d\phi = 2 / (1 + r^2) dr$
$\sin(\phi)^{n - 1} dL \, d\phi = 2 \sin(\phi)^{n - 1} dL / (1 + r^2) dr$
$dA_{\text{sphere}} = 2 \sin(\phi)^{n - 1} r^{-n + 1} / (1 + r^2) dA_{\text{plane}}$
$\sin(\phi) = \sin(2 \tan^{-1}(r)) = 2r / (1 + r^2)$
$\therefore$ The PDF of the distribution on the plane should be $p(x) = c \, (2r / (1 + r^2))^{n - 1} r^{-n + 1} / (1 + r^2)$ for some normalization factor $c$.
This is when I realize that I can't easily integrate the PDF and find the inverse function for the CDF, which I need to do to find $f$. Mathematica couldn't find a closed form formula for most values of $n$. Maybe there a pattern to the integrals and a quick way to numerically evaluate the inverse function?