Skip to main content
3 of 9
added 1447 characters in body
Thomas Steinke
  • 2.8k
  • 20
  • 29

My understanding of your question is that you want to sample a point uniformly (i.e. Lebesgue measure) from the ball $$B(d,r)=\left\{ x \in \mathbb{R}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ I will give two algorithms for this. For simplicity, I assume $r=1$, as you can always scale up.


Algorithm 1: Rejection Sampling
  1. Sample $U_1, \cdots, U_d \in [-1,1]$ independently and uniformly at random.
  2. Let $x = (U_1, \cdots, U_d)$.
  3. If $\|x\|_2 \leq 1$, return $x$; else GOTO 1.
This algorithm is simple. We just use the fact that it is easy to sample from the "box" $[-1,1]^d$ that contains the unit ball. And, by rejecting samples that fall outside the ball, we are able to sample from the desired conditional distribution. The downside of this algorithm is that it takes time exponential in $d$, as the probability of a random point from the box being in the ball is decaying exponentially with $d$.
Algorithm 2: Scaled Spherical Sample
  1. Sample $G_1, \cdots, G_d \in \mathbb{R}$ independently from a standard Gaussian distribution ($\mathcal{N}(0,1)$).
  2. Let $x = (G_1, \cdots, G_d)$ and $y=x/\|x\|_2$.
  3. Sample $U \in [0,1]$ uniformly at random.
  4. Return $z=y \cdot U^{1/d}$.
Since the Gaussian distribution is spherically symmetric, $y$ is a uniformly random unit vector. Now we scale down $y$ by $U^{1/d}$ so that, instead of uniformly sampling from the sphere, we uniformly sample from the ball. To show that this is the right scaling, we just need to show that $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}}.$$ We have $$\frac{\text{Volume of $d$-ball of radius $t$}}{\text{Volume of $d$-ball of radius $1$}} = \frac{\frac{\pi^{d/2}}{\Gamma(d/2+1)} t^{d}}{\frac{\pi^{d/2}}{\Gamma(d/2+1)} 1^{d}} = t^{d}$$ and $$\mathbb{P}\left[\|z\|_2 \leq t\right] = \mathbb{P}\left[U^{1/d} \leq t\right] = \mathbb{P}\left[U \leq t^d\right] = t^d,$$ as required. Finally note that we can efficiently sample standard Gaussians using the Box-Muller transform.

EDIT: I see from the comments that you probably want an integer point -- that is, uniformly random from $$B_\mathbb{Z}(d,r)\left\{ x \in \mathbb{Z}^d : \|x\|_2 = \sqrt{\sum_{i=1}^d x_i^2} \leq r \right\}.$$ Peter Shor's answer covers this case. Below is an alternative approach.

The approach is to sample one coordinate at a time. Define $f(d,r^2)=\left|B_\mathbb{Z}(d,r)\right|$. We first use $f$ to sample from $B_\mathbb{Z}(d,r)$ and then show how to compute $f$.

We can partition $$B_\mathbb{Z}(d,r) = \bigcup_{a \in \mathbb{Z}} \{a\} \times B_\mathbb{Z}\left(d-1,\sqrt{r-a^2}\right).$$ Thus we can (i) sample $a \in \mathbb{Z}$ from the appropriate marginal distribution and (ii) recursively sample from $B_\mathbb{Z}\left(d-1,\sqrt{r-a^2}\right)$. Combining the two gives a uniform sample from $B_\mathbb{Z}(d,r)$.

To sample $a$, we compute $f(d,r^2-a^2)$ for all $a \in [-r,r] \cap \mathbb{Z}$ and then scale these values to gain a probability mass function for $a$.

All that remains is to determine how to compute $f$. We can use the partition to obtain a recursive formula $$f(d,r^2) = \sum_{a \in [-r,r] \cap \mathbb{Z}} f(d-1,r^2-a^2).$$ The base case of the recursion is $f(0,r^2)=1$ if $r\geq0$ and $f(d,r)=0$ if $r<0$. So we can compute all necessary values of $f$ using dynamic programming with a lookup table of size $O(r^2d)$.

Packing this into a single algorithm:

Algorithm 3: Recursion
  1. Allocate a lookup table for $f : \{0, \cdots, d\} \times \{0, \cdots, r^2\} \to \mathbb{Z}$.
  2. Initialize $f(0,s)=1$ for $s = 0 \cdots r^2$.
  3. For $m = 1 \cdots d$:
  4.    For $s = 0 \cdots r^2$:
  5.       Set $f(m,s) = \sum_{a = \lceil - \sqrt{s} \rceil}^{\lfloor \sqrt{s} \rfloor} f(m-1,\lfloor s-a^2 \rfloor).$
  6. Set $s_1=r^2$ and allocate space for $x \in \{-r,r\}^d$.
  7. For $i = 1 \cdots d$:
  8.    Sample $U_i \in \{1, \cdots, f(d+1-i,s_i)\}$ uniformly at random.
  9.    Pick the largest $x_i \in \{\lceil - \sqrt{s_i} \rceil, \cdots, \lfloor \sqrt{s_i} \rfloor\}$ such that $$\sum_{a = \lceil - \sqrt{s_i} \rceil}^{x_i} f(d-i,\lfloor s_i-a^2 \rfloor) \leq U_i.$$
  10.    Set $s_{i+1}=s_i-x_i^2$.
  11. Return $x$.
Thomas Steinke
  • 2.8k
  • 20
  • 29