How to generate random points in $\ell_p$ balls?

Question

How do I feasibly generate a random sample from an $n$-dimensional $\ell_p$ ball? Specifically, I'm interested in $p=1$ and large $n$. I'm looking for descriptions analogous to the statement for $p=2$: Take $n$ standard gaussian random variables and normalize.

Answer 1

For arbitrary p, this paper does exactly what you want. Specifically, pick $X_1,\ldots,X_n$ independently with density proportional to $\exp(-|x|^p)$, and $Y$ an independent exponential random variable with mean 1. Then the random vector $$\frac{(X_1,\ldots,X_n)}{(Y+\sum |X_i|^p)^{1/p}}$$ is uniformly distributed in the unit ball of $\ell_p^n$.

The paper also shows how to generate certain other distributions on the $\ell_p^n$ ball by modifying the distribution of $Y$.

Answer 2

---

I'll assume that you're looking for a uniformly chosen random point in the ball, since you didn't state otherwise. For p=1, you're asking for a uniform random point in the cross polytope in n dimensions. That is the set

$ C_n = \{ x_1, x_2, \ldots, x_n \in \mathbb{R} : |x_1| + \cdots + |x_n| \le 1 \}. $

By symmetry, it suffices to pick a random point $(X_1, \ldots, X_n)$ from the simplex

$ S_n = \{ x_1, x_2, \ldots, x_n \in \mathbb{R}^+ : x_1 + \cdots + x_n \le 1 \}$

and then flip $n$ independent coins to attach signs to the $x_i$.

From Devroye's book Non-uniform random variable generation (freely available on the web at the link above, see p. 207 near the beginning of Chapter 5), we can pick a point in the simplex uniformly at random by the following procedure:

• let $U_1, \ldots, U_n$ be iid uniform(0,1) random variables
• let $V_1, \ldots, V_n$ be the $U_i$ reordered so that $V_1 \le V_2 \le \cdots \le V_n$ (the "order statistics"); let $V_0 = 0, V_{n+1} = 1$
• let $X_i = V_i - V_{i-1}$

So do this to pick the absolute values of the coordinates of your points; attach signs chosen uniformly at random, and you're done.

This of course relies on the special structure of balls in $\ell^1$; I don't know how to generalize it to arbitrary $p$.

Answer 3

[Here is a draft I'm working on and sent to a journal. Posting it here I am guessing violates current publication practices, but would love to get feedback on this here]

Draft: https://drive.google.com/drive/folders/1W7y_Nc-4XjuI-pfkTdCXY6QATMnzY3YX?usp=sharing

---

Let ${ p \in [1, \infty) }.$ Now ${ \lVert \ldots \rVert _{p} }$ is a norm on ${ \mathbb{R} ^n },$ and we write ${ \ell _{p} ^{n} }$ for the space ${ \ell _{p} ^{n} := (\mathbb{R} ^n, \lVert \ldots \rVert _p) . }$

The problem is to find a random vector uniform over unit ball ${ B _p := \lbrace x \in \mathbb{R} ^n : \lVert x \rVert \leq 1 \rbrace. }$

Consider first the positive part of unit ball ${ B _p ^{+} := \lbrace x \in \mathbb{R} ^n : \text{each } x _i \geq 0, \lVert x \rVert _{p} \leq 1 \rbrace. }$

Any point in ${ B _{p} ^+ }$ can be written as $${ \left( \frac{x _1}{(\sum _{1} ^{n} x _i ^{p} + x _{n+1}) ^{\frac{1}{p}}} , \ldots, \frac{x _n}{(\sum _{1} ^{n} x _i ^{p} + x _{n+1}) ^{\frac{1}{p}} } \right) }$$ with each ${ x _i \geq 0 }.$
Intuitively any vector ${ y = \left( \frac{x _1}{(\sum _{1} ^{n} x _i ^{p} + x _{n+1}) ^{\frac{1}{p}}} , \ldots, \frac{x _n}{(\sum _{1} ^{n} x _i ^{p} + x _{n+1}) ^{\frac{1}{p}} } \right) }$ with each ${ x _i \geq 0 }$ is parallel to ${ (x _1, \ldots, x _{n}) }$ with $p-$norm ${ \lVert y \rVert _p }$ ${ = \left( \frac{\sum _{1} ^{n} x _{i} ^{p}}{\sum _{1} ^{n} x _{i} ^{p} + x _{n+1} } \right) ^{\frac{1}{p}} \leq 1 },$ and if ${ (x _1, \ldots, x _n) \neq 0 }$ is fixed the $p-$norm ${ \lVert y \rVert _{p} }$ can be set to take any value in ${ (0,1] }$ by appropriate choice of ${ x _{n+1} \geq 0 }.$

Say ${ X _1, X _2, \ldots \geq 0 }$ are independent positive random variables with densities ${ f _1, f _2, \ldots }$ respectively. We will pick "nice" ${ f _i }$s such that the random vector $${ \left( \frac{X _1}{(\sum _{1} ^{n} X _i ^{p} + X _{n+1}) ^{\frac{1}{p}}} , \ldots, \frac{X _n}{(\sum _{1} ^{n} X _i ^{p} + X _{n+1}) ^{\frac{1}{p}} } \right) }$$ is uniform over ${ B _{p} ^{+} }.$

Since change-of-density theorem doesnt apply to the mapping ${ (0, \infty) ^{n+1} \to \text{int}(B _{p} ^{+}), }$ ${ (x _1, \ldots, x _{n+1}) \mapsto \left( \frac{x _1}{(\sum _{1} ^{n} x _i ^{p} + x _{n+1}) ^{\frac{1}{p}}} , \ldots, \frac{x _n}{(\sum _{1} ^{n} x _i ^{p} + x _{n+1}) ^{\frac{1}{p}} } \right) }$ due to reduction in dimension, we can instead look at $${ F : (0,\infty) ^{n+1} \to \text{int}(B _{p} ^{+}) \times (0, \infty), }$$ $${ (x _1, \ldots, x _{n+1}) \mapsto (y _1, \ldots, y _{n+1}) = \left( \frac{x _1}{(\sum _{1} ^{n} x _i ^{p} + x _{n+1}) ^{\frac{1}{p}}} , \ldots, \frac{x _n}{(\sum _{1} ^{n} x _i ^{p} + x _{n+1}) ^{\frac{1}{p}} } , (\sum _{1} ^{n} x _i ^{p} + x _{n+1}) ^{\frac{1}{p}} \right) }$$ which can be inverted as $${ F ^{-1} : (y _1, \ldots, y _{n+1}) \mapsto (x _1, \ldots, x _{n+1}) = (y _1 y _{n+1}, \ldots, y _{n} y _{n+1}, (y _{n+1})^{p} - (y _1 y _{n+1})^{p} - \ldots - (y _n y _{n+1})^{p} ) }.$$

Now random vector ${ (Y _1, \ldots, Y _{n+1}) = F(X _1, \ldots, X _{n+1}) \in \text{int}(B _{p} ^{+}) \times (0, \infty) }$ has density $${ f _{Y} (y) = f _{X} (F ^{-1} (y)) \vert \det D(F ^{-1}) _{y} \vert . }$$
The derivative $${ D(F ^{-1}) _{y} = \begin{pmatrix} y _{n+1} &0 &\ldots &0 &y _{1} \\ 0 &y _{n+1} &\ldots &0 &y _{2} \\ \vdots &\vdots &\ddots &\vdots &\vdots \\ 0 &0 &\ldots &y _{n+1} &y _{n} \\ \partial _{1} x _{n+1} &\partial _{2} x _{n+1} &\ldots &\partial _{n} x _{n+1} &\partial _{n+1} x _{n+1} \end{pmatrix}, }$$ where entries $${ \partial _{i} x _{n+1} = (y _{n+1}) ^{p} (-p y _{i} ^{p-1}) \text{ for } i = 1, \ldots, n}$$ and $${ \partial _{n+1} x _{n+1} = p (y _{n+1}) ^{p-1} (1- (y _1) ^{p} - \ldots - (y _{n}) ^{p}) }.$$ So successively adding row ${ i }$ times ${p (y _{n+1}) ^{p-1} (y _i) ^{p-1} }$ to last row, for ${ i = 1, \ldots, n },$ gives its determinant to be \begin{align*} \det D(F ^{-1}) _{y} &= (y _{n+1}) ^{n} [\partial _{n+1} x _{n+1} + \sum _{1} ^{n} y _i (p (y _{n+1}) ^{p-1} (y _i) ^{p-1} )] \\ &= p (y _{n+1}) ^{n+p-1}. \end{align*}

So for ${ y \in \text{int}(B _{p} ^{+}) \times (0, \infty) }$ density $${ f _{Y} (y) = p(y _{n+1}) ^{n+p-1} f _{1} (y _{1} y _{n+1}) \ldots f _{n} (y _{n} y _{n+1}) f _{n+1} ((y _{n+1})^{p}(1- y _{1} ^{p} - \ldots - y _{n} ^{p} )) . }$$

Hence the vector of interest ${ (Y _1, \ldots, Y _{n}) \in \text{int}(B _{p} ^{+}) }$ has density \begin{equation}\label{eq:changeofdensity}\begin{split} &\quad f _{Y _1, \ldots, Y _{n}} (y _1, \ldots, y _n) \\ &= \int _{0} ^{\infty} f _{Y _1, \ldots, Y _{n+1}} (y _1, \ldots, y _{n+1}) \, dy _{n+1} \\ &= p \int _{0} ^{\infty} z ^{n+p-1} f _{1} ({\color{red}{y _{1}}} z) \ldots f _{n} ({\color{red}{y _{n}}} z) f _{n+1} (z ^{p}(1-{\color{red}{y _{1} ^{p}}} - \ldots - {\color{red}{y _{n} ^{p}}})) \, dz . \end{split}\end{equation} We want ${ f _{i} }$s to be such that above density is constant over ${ ({\color{red}{y _1}}, \ldots, {\color{red}{y _{n}}}) \in \text{int}(B _{p} ^{+}) }.$

This will happen if $${ f _{1} ({\color{red}{y _{1}}} z) \ldots f _{n} ({\color{red}{y _{n}}} z) f _{n+1} (z ^{p}(1-{\color{red}{y _{1} ^{p}}} - \ldots - {\color{red}{y _{n} ^{p}}})) }$$ is a function of ${ z }$ alone.

If densities ${ f _{1}(t) , \ldots, f _{n}(t) }$ are all proportional to ${ e ^{-t ^p} }$ (i.e. ${ f _{1} (t) = \ldots = f _{n} (t) = \frac{1}{\Gamma(1+\frac{1}{p})} e ^{-t ^{p}} }$ for ${ t \geq 0 }$) and density ${ f _{n+1} (t) }$ is proportional to ${ e ^{-t} }$ (i.e. ${ f _{n+1} (t) = e ^{-t} }$ for ${ t \geq 0 }$), then $${ f _{1} ({\color{red}{y _{1}}} z) \ldots f _{n} ({\color{red}{y _{n}}} z) f _{n+1} (z ^{p}(1-{\color{red}{y _{1} ^{p}}} - \ldots - {\color{red}{y _{n} ^{p}}})) = \frac{1}{\Gamma(1+\frac{1}{p}) ^{n} } e ^{-z ^{p}}, }$$ as needed. This gives:

Thm [Sampling uniformly from ${ B _p ^{+} }$]: Let ${ p \in [1, \infty) }.$ If ${ X _1, \ldots, X _n, X _{n+1} }$ are independent positive random variables, with densities proportional to ${ e ^{-t ^p}, \ldots, e ^{-t ^p} , e ^{-t} }$ (for ${ t \geq 0 }$) respectively, the random vector $${ (Y _1, \ldots, Y _n) = \left( \frac{X _1}{(\sum _{1} ^{n} X _i ^{p} + X _{n+1}) ^{\frac{1}{p}}} , \ldots, \frac{X _n}{(\sum _{1} ^{n} X _i ^{p} + X _{n+1}) ^{\frac{1}{p}} } \right) }$$ is uniform over ${ B _{p} ^{+} = \lbrace x \in \mathbb{R} ^{n} : \text{each } x _i \geq 0, \lVert x \rVert _{p} \leq 1 \rbrace }.$

If a random vector ${ (Y _1, \ldots, Y _n) }$ is uniform over ${ B _{p} ^{+} },$ the vector ${ ((-1) ^{Z _1} Y _1, \ldots, (-1) ^{Z _n} Y _n) }$ where ${ Z _1, \ldots, Z _n }$ are independent ${ \text{Unif}\lbrace 0,1 \rbrace }$ variables is uniform over the ball ${ B _p = \lbrace x \in \mathbb{R} ^{n} : \lVert x \rVert _{p} \leq 1 \rbrace }.$ In the context of above theorem, $${ ((-1) ^{Z _1} Y _1, \ldots, (-1) ^{Z _n} Y _n) = \left( \frac{(-1) ^{Z _1} X _1}{(\sum _{1} ^{n} \vert (-1) ^{Z _i} X _i \vert ^{p} + X _{n+1} ) ^{\frac{1}{p}}} , \ldots, \frac{(-1) ^{Z _n} X _n}{(\sum _{1} ^{n} \vert (-1) ^{Z _i} X _i \vert ^{p} + X _{n+1} ) ^{\frac{1}{p}} } \right), }$$ and the variables ${ (-1) ^{Z _1} X _1, \ldots, (-1) ^{Z _n} X _n , X _{n+1} }$ are independent with densities proportional to ${ e ^{-\vert t \vert ^{p}} , \ldots, e ^{-\vert t \vert ^{p}} }$ (for ${ t \in \mathbb{R} }$) and ${ e ^{-t} }$ (for ${ t \geq 0 }$) respectively. This gives, as in the paper https://arxiv.org/abs/math/0503650 linked above due to Barthe, Guedon, Mendelson and Naor:

Thm [Sampling uniformly from ${ B _p }$]: Let ${ p \in [1, \infty) }.$ If ${ X _1, \ldots, X _{n+1} }$ are independent random variables, with ${ X _1, \ldots, X _{n} }$ having densities proportional to ${ e ^{- \vert t \vert ^{p}} }$ (for ${ t \in \mathbb{R} }$) and ${ X _{n+1} }$ having density ${ e ^{-t} }$ (for ${ t \geq 0 }$), the random vector $${ (Y _1, \ldots, Y _n) = \left( \frac{X _1}{(\sum _{1} ^{n} \vert X _i \vert ^{p} + X _{n+1}) ^{\frac{1}{p}}} , \ldots, \frac{X _n}{(\sum _{1} ^{n} \vert X _i \vert ^{p} + X _{n+1}) ^{\frac{1}{p}} } \right) }$$ is uniform over ${ B _p = \lbrace x \in \mathbb{R} ^{n} : \lVert x \rVert _{p} \leq 1 \rbrace }.$

---

Remark: The same method $${ f _{Y _1, \ldots, Y _{n}} (y _1, \ldots, y _n) = p \int _{0} ^{\infty} z ^{n+p-1} f _{1} ({\color{red}{y _{1}}} z) \ldots f _{n} ({\color{red}{y _{n}}} z) f _{n+1} (z ^{p}(1-{\color{red}{y _{1} ^{p}}} - \ldots - {\color{red}{y _{n} ^{p}}})) \, dz }$$ but different choices of densities ${ f _i }$s gives different densities on ${ B _p ^+ }$ and ${ B _p , }$ as outlined in the draft.