11
$\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 '09 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 '09 at 17:38
24
$\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
  • 3
    $\begingroup$ This is the natural answer to the problem. It directly generalizes the Gaussian point method for $p=2$. $\endgroup$ Dec 17 '09 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 '18 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 '18 at 21:38
10
$\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 '09 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 '09 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 '09 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 '09 at 17:50
  • $\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 '20 at 14:37

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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