Projecting a distribution onto the unit sphere in an arbitrary norm

So we represent the Dirichlet distribution as the projection of the $$d$$ independent gammas (on $$R_+^d$$ onto the unit simplex, and we arrive at that through the $$L_1$$ norm. That is, divide $${\bf x} \in R_+^{d}$$ by $$r = \sum x_i$$ to push data onto the unit simplex, $$S_{1}^{d-1}$$. Then we integrate out $$r$$, and are left with Dirichlet (or generalized dirichlet, if using variable rate gammas).

I'm trying to set up the Jacobian of the transformation between $${\bf x} \in R_+^d$$ and $${\bf y} \in S_{\omega}^{d-1}$$, where $$\omega$$ specifies some arbitrary norm. If we want to narrow it down, then in the family of $$p$$-norms. That is, $$r = \left(\sum x_i^p\right)^{1/p}$$ So you can imagine, $$p = 1$$ defines the simplex, $$p = 2$$ defines the hypersphere. My problem is specifically defined on the infinity norm; $$\lim_{p\to\infty} (\sum x_i^p)^{1/p} = \max x_i$$

The determinant of the Jacobian should still be $$r^{d-1}$$, but I'm not certain of that. It's really the ordering of variables that's messing me up; I need to end up with a random vector $$(r,y_{(1)},\ldots,y_{(d-1)})$$, and which dimension falls out is the max, but that's random.

Is it sufficient to show the determinant of the jacobian is the same no matter which dimension falls out?

EDIT: To be clear, the standard transformation of the is something like: \begin{aligned} T(x_1,\ldots,x_d) &= (r, y_1,\ldots,y_{d-1})\\ \implies T^{-1}(r,{\bf y}) &= \left(ry_1,\ldots,ry_{d-1},r(1-\sum_{i=1}^{d-1}y_i^p)^{1/p}\right) \end{aligned} But when I take the determinant of the Jacobian of that transformation, I end up with something like $$r^{d-1}\left[(1 - \sum_{i = 1}^{d-1}y_i^p)^{1/p} + \sum_{i = 1}^{d-1}y_i^p / (1 - \sum_{l = 1}^{d-1}y_l^p)^{1/p - 1}\right].$$ If I take $$p\to\infty$$, if anything but $$y_d = 1$$, then the denominator goes to 0.

• Your definition of $T^{-1}$ doesn't make sense, because the right hand side depends on the $x_i.$ That might be what's getting in your way. I suspect you will need to compute the Jacobian in a different manner for $p=\infty$ (since it's not everywhere differentiable), but that's relatively easy because locally the inverse transformation is multilinear.
– whuber
Commented Apr 23, 2021 at 19:47

Suppose $$(X_1,\ldots,X_d)$$ is a positive continuous random variable. Write its joint probability element as

$$f(x_1,\ldots,x_d)\,\mathrm{d}x_1 \cdots \mathrm{d}x_{d-1} \mathrm{d}x_d$$

so that the density function is $$f(x_1,\ldots, x_d).$$

Leaving the case of the $$L^\infty$$ norm for later (it is special), consider $$0\le p\lt \infty.$$ You ask about the density of the random variable

$$(Y_1,\ldots,Y_{d-1},R) = (X_1/R, \ldots, X_{d-1}/R, (X_1^p+\cdots+X_d^p)^{1/p}).$$

This is an invertible transformation because $$X_i=R\,Y_i$$ is defined for $$i=1,2,\ldots,d-1$$ and $$X_d = R\left(1 - Y_1^p-\cdots-Y_{d-1}^p\right)^{1/p}$$ also is unique.

Emulating the derivation of the Dirichlet distribution using the algebra of differential forms, compute

$$\mathrm{d}x_i = r\,\mathrm{d}y_i + y_i\,\mathrm{d}r,\ i=1,2,\ldots, d-1$$

and

\begin{aligned} px_d^{p-1}\,\mathrm{d}x_d &= d\left(x_d^p\right) = \mathrm d\left(r^p-x_1^p-\cdots-x_{d-1}^p\right)\\ &= p\left(r^{p-1}\,\mathrm{d}r - x_1^{p-1}\,\mathrm{d}x_1 - \cdots - x_{d-1}^{p-1}\,\mathrm{d}x_{d-1}\right). \end{aligned}

Plug these in to compute the old differential in terms of the new variables and their differentials:

\begin{aligned} \mathrm{d}x_1\wedge \cdots \wedge \mathrm{d}x_d &= \left(px_d^{p-1}\right)^{-1}\mathrm{d}x_1\wedge \cdots \wedge \mathrm{d}x_{d-1}\wedge \mathrm{d}(x_d^p)\\ &= \left(x_d^{p-1}\right)^{-1}\mathrm{d}x_1\wedge \cdots \wedge \mathrm{d}x_{d-1}\wedge \left(r^{p-1}\,\mathrm{d}r - x_1^{p-1}\,\mathrm{d}x_1 - \cdots - x_{d-1}^{p-1}\,\mathrm{d}x_{d-1}\right)\\ &= \left(x_d^{p-1}\right)^{-1}\mathrm{d}x_1\wedge \cdots \wedge \mathrm{d}x_{d-1}\wedge r^{p-1}\,\mathrm{d}r \\ &= \left(x_d^{p-1}\right)^{-1} \left(r\,\mathrm{d}y_1 + y_1\,\mathrm{d}r\right) \wedge \cdots \wedge \left(r\,\mathrm{d}y_{d-1} + y_{d-1}\,\mathrm{d}r\right) \wedge r^{p-1}\,\mathrm{d}r\\ &= r^{d-1}\left(\frac{r^{p-1}}{x_d^{p-1}}\right)\,\mathrm{d}y_1\wedge \cdots \wedge \mathrm{d}y_{d-1} \wedge \mathrm{d}r\\ &= r^{d-1}\left(1 - y_1^p - \cdots - y_{d-1}^p\right)^{1/p-1}\,\mathrm{d}y_1\wedge \cdots \wedge \mathrm{d}y_{d-1} \wedge \mathrm{d}r. \end{aligned}

The coefficient of the differential form is the Jacobian and (therefore) the density of $$(Y_1,\ldots, Y_{d-1})$$ is

\begin{aligned} &g(y_1,\ldots, y_{d-1}) \\ &= \int_0^\infty f\left(ry_1, \ldots, ry_{d-1}, r\left(1-y_1^p-\cdots-y_{d-1}^p\right)^{1/p}\right)\, r^{d-1}\left(1 - y_1^p - \cdots - y_{d-1}^p\right)^{1/p-1}\,\mathrm{d}r. \end{aligned}

The difficulty when $$p=\infty,$$ where $$r=\max(x_i),$$ is that $$(y_1,\ldots, y_{d-1}) = (x_1/r,\ldots, x_{d-1}/r)$$ are no longer coordinates for the image of the projection. These coordinates work only in the subset where $$x_d$$ is the largest.

A simple approach is to decompose the distribution according to the events determined by which of the $$X_i$$ is the largest. Because the joint distribution is continuous, the chance of a tie for largest is zero, so may neglect that possibility. Let $$K$$ be the random variable given by the index of the largest of the $$X_i.$$

Conditioning first on $$K=d,$$ we have $$R=\max(X_i)=X_d$$ and $$Y_i = X_i/R = X_i/X_d$$ for $$i=1,2,\ldots, d-1.$$ Computing as before we (easily) obtain

$$\mathrm{d}x_1\wedge \cdots \wedge \mathrm{d}x_{d} = r^{d-1}\, \mathrm{d}y_1\wedge \cdots \wedge \mathrm{d}y_{d-1}\wedge \mathrm{d}r.$$

Comparing this to the $$L^p$$ solution obtained previously, we may offer an intuitive (hand-waving) argument for this formula (as suggested in the question): as $$p$$ grows large, the factor $$(1-y_1^p-\cdots -y_{d-1}^p)^{1/p-1}$$ approaches $$x_d^{-1} = r^{-1},$$ giving $$r^d(r^{-1})= r^{d-1}$$ for the Jacobian.

This analysis applies with almost no change when conditioning on other values of $$K:$$ you only have to take the proper sign when integrating over $$r$$ in each case to assure the result is a positive density. In every case the Jacobian still equals $$r^{d-1}.$$ The full distribution of the projection is an equal mixture of the $$d$$ conditional distributions. Each of the $$d$$ components is supported on a different face of the unit cube $$I_d = (0,1]^d\subset \mathbb{R}^d;$$ the component corresponding to $$K=k$$ is supported on $$\partial_k^{+} I_d = \{(y_1,\ldots,y_d)\mid 0\lt y_i\le 1;\ y_k=1\} \subset I_d.$$