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.


Bin Fu

  • $\begingroup$ What do you mean by "grid point"? Do you just want to sample uniformly from the L2 ball? $\endgroup$ – Thomas Feb 7 '16 at 6:07
  • $\begingroup$ If you do not require lattice points, then you can do the following (for sampling uniformly from the L2 ball): Pick a random direction with any spherically symmetric distribution, e.g. picking d coordinates using identical independent normal distributions and using the resulting direction, and then choose the radius, which needs to be chosen with a density for radius $r$ that is proportional to $r^{d-1}$. I suppose you are actually are interested in lattice points (integer points) that lie within the hyperball, but as Thomas requested, please clarify. $\endgroup$ – Rahul Savani Feb 7 '16 at 8:20
  • 1
    $\begingroup$ What distribution? Uniform? $\endgroup$ – usul Feb 7 '16 at 14:21
  • $\begingroup$ Yes, I am interested in lattice point as you pointed out. Your algorithm makes some sense. We need to control time and bias. $\endgroup$ – Bin Fu Feb 7 '16 at 14:51
  • $\begingroup$ I look for random lattice point under uniform distribution, but allow a small bias. Currently I put the ball in a cube and generate random lattice points until one of them is in the ball. It still takes exponential time. $\endgroup$ – Bin Fu Feb 7 '16 at 15:04

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.


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.

  • $\begingroup$ I like Peter and your suggestions. There may be a little bias that I will check it $\endgroup$ – Bin Fu Feb 7 '16 at 19:48
  • $\begingroup$ Unless I made a mistake, the algorithms I suggested should give perfectly uniform outputs. Peter's first algorithm also outputs a perfectly uniform point, but his second one only outputs an approximately uniform point. $\endgroup$ – Thomas Feb 7 '16 at 19:52
  • $\begingroup$ It seems that function f(.) depends on the center if the original ball center $\endgroup$ – Bin Fu Feb 7 '16 at 21:20
  • $\begingroup$ is not at a lattice point. Otherwise, it is uniform. $\endgroup$ – Bin Fu Feb 7 '16 at 21:21
  • $\begingroup$ I am assuming the ball is centered at 0. Is this not the case? $\endgroup$ – Thomas Feb 7 '16 at 21:23

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.

  • $\begingroup$ This is an exponential time algorithm with dimension number d. I am extremely interested in the case when d is large and hope see the $\endgroup$ – Bin Fu Feb 7 '16 at 6:55
  • $\begingroup$ Want to see a polynomial time solution. $\endgroup$ – Bin Fu Feb 7 '16 at 6:56

Your Answer

By clicking "Post Your Answer", you acknowledge that you have read our updated terms of service, privacy policy and cookie policy, and that your continued use of the website is subject to these policies.

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