14
$\begingroup$

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.

$\endgroup$
2
  • $\begingroup$ There's an obvious map from the $\ell_2$ ball to the $\ell_1$ ball. Does that not work for you? What kind of distribution do you want? $\endgroup$ Dec 17, 2009 at 17:27
  • $\begingroup$ I'm interested in the uniform distribution. If by "obvious map" you mean rescaling, then that wouldn't work. $\endgroup$
    – Mitch
    Dec 17, 2009 at 17:38

3 Answers 3

31
$\begingroup$

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$.

$\endgroup$
3
  • 5
    $\begingroup$ This is the natural answer to the problem. It directly generalizes the Gaussian point method for $p=2$. $\endgroup$ Dec 17, 2009 at 19:40
  • $\begingroup$ Is there a way to sample exactly from $\exp(-|x|^p)$ for $p>2$? or should I just use slice sampling or other MCMC method? $\endgroup$
    – gwding
    Feb 1, 2018 at 20:30
  • $\begingroup$ @gwding: I'm not an expert on random number generation, and I'm really not sure what the done thing is when (as in this case) there's no simple formula for the cdf. $\endgroup$ Feb 1, 2018 at 21:38
12
$\begingroup$

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$.

$\endgroup$
5
  • $\begingroup$ C_n and S_N are of course sets; I can't figure out how to get the braces to show up. $\endgroup$ Dec 17, 2009 at 17:37
  • $\begingroup$ Thanks for the reference. For $p=1$ this was what I thought, but I didn't have a reference. I'm still interested in p>2 though. $\endgroup$
    – Mitch
    Dec 17, 2009 at 17:45
  • $\begingroup$ Fixed the braces. If you're curious how I did this, I'll put up a note on tea.mathoverflow.net $\endgroup$ Dec 17, 2009 at 17:48
  • $\begingroup$ David, thanks for fixing the notation. Mitch, I'd love an answer for other p as well. It seems like a problem someone should have considered but quick searching isn't finding any answers other than the obvious one of picking a random point in the unit cube and throwing it out if it doesn't work, which for large n gets very bad. $\endgroup$ Dec 17, 2009 at 17:50
  • 1
    $\begingroup$ With this process the L1 norm of X is always $V_{n+1} - V_0 = 1$ so this simulates a random point on the $l_1$ sphere, not the ball $\endgroup$
    – Gericault
    May 18, 2020 at 14:37
0
$\begingroup$

[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.

$\endgroup$
3
  • 1
    $\begingroup$ It doesnt violate publication practices, but it's not good form on MO to just link to a pdf from an answer with no explanation. If this is intended to answer the question, please outline what your document says. If this posting is only intended to attract feedback on a draft then that's really not what MO is for, I'm afraid. $\endgroup$
    – David Roberts
    Oct 15, 2023 at 1:46
  • 1
    $\begingroup$ Thanks, I will try to outline relevant part of the argument here. $\endgroup$ Oct 15, 2023 at 5:37
  • 1
    $\begingroup$ Thanks, that's much more informative! $\endgroup$
    – David Roberts
    Oct 15, 2023 at 9:32

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.