Random grid point in a d-dimensional ball

Question

I would like to know if there is any standard algorithm to generate a random grid point inside a d-dimensional ball with a given radius r.

Thanks

Bin Fu

Answer 1

Let $B$ be the ball you want to sample in. Let $B'$ be the ball centered at the same point, but larger in radius by $\sqrt{d}$. Sample a random point in $B'$ Round off to the nearest integer point. Now, if this integer point is in $B$, accept it. Otherwise, discard it and try again.

This works because every point of the original ball is contained in a unit cube in ball $B'$.

When will this be efficient? You want the ball of $r+ \sqrt{d}$ to be not too much larger than the ball of radius $r$. The ratio of these distances is

$$ \left(\frac{r+\sqrt{d}}{r}\right)^d = \left(1 + \frac{\sqrt{d}}{r}\right)^d \leq \exp\left({\frac{d^{3/2}}{r}}\right).$$

So it works well if $r > d^{3/2}$.

How do you find a random point in $B'$? Thomas's answer give an algorithm for this.

---

If $r < d^{3/2}$, use the Metropolis algorithm. Start with a point $P$ in the ball—the lattice point closest to the center is easiest. Now, repeatedly perform the following steps.

1. Choose a random coordinate.

2. Keeping the remaining coordinates fixed, compute the lowest and highest integer values for this coordinate that are still in the ball.

3. Pick a random value $v$ for this coordinate, and set this coordinate to $v$ in $P$.

This is guaranteed to converge, and there are theorems that will tell you how quickly it will converge ... I am fairly sure that this case is one that converges in polynomial time.

Answer 2

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.

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

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

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

All that remains is to determine how to compute $f$. We can use the above 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$. 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

• Input: $d,r^2 \in \mathbb{N}$.
• Output: $x \in B_\mathbb{Z}(d,r)$ uniformly random.

1. Allocate a lookup table for $f : \{0, \cdots, d\} \times \{0, \cdots, r^2\} \to \mathbb{N}$.
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, s-a^2 ).$
6. Set $s_1=r^2$ and allocate space for $x \in \{\lceil -r \rceil, \cdots ,\lfloor r \rfloor\}^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, s_i-a^2 ) \leq U_i.$$
10.    Set $s_{i+1}=s_i-x_i^2$.
11. Return $x \in B_\mathbb{Z}(d,r)$.

The running time of this algorithm is $O(r^3d)$. Assuming the $U_i$s are uniform and independent, it outputs a uniformly random point from $B_\mathbb{Z}(d,r)$.

---

Before EDIT: 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.

Answer 3

Just off the top of my head if you want to do it repeatedly and can afford the memory: 1) Enumerate your "grid points" in advance with each possible coordinate d-tuple assigned to a number from $0$ to $n-1$, where there are $n$ grid points. 2) Pick a random number from $0$ to $n-1$ and look up which coordinate d-tuple it corresponds to. This would use a bit of memory but make each random sampling be as fast as generating a random integer and reading that index in an array.