<p>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.</p>

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

<p>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 </p>

<p>$$ \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).$$</p>

<p>So it works well if $r &gt; d^{3/2}$.</p>

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

<hr>

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

<ol>
<li><p>Choose a random coordinate.</p></li>
<li><p>Keeping the remaining coordinates fixed, compute the lowest and highest integer values for this coordinate that are still in the ball.</p></li>
<li><p>Pick a random value $v$ for this coordinate, and set this coordinate to $v$ in $P$.</p></li>
</ol>

<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.</p>

CC BY-SA 3.0
2016-02-07T15:05:33.157Z
2016-02-07T18:46:50.520Z
2016-02-07T18:46:50.520Z
Peter Shor 
<p><b>EDIT:</b> 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\}.$$
<a href="https://cstheory.stackexchange.com/a/33755/2022">Peter Shor's answer</a> covers this case. Below is an alternative approach.</p>

<p>The approach is to sample one coordinate at a time. </p>

<p>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)$.</p>

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

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

<p>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)$.</p>

<p>Packing this into a single algorithm:</p>

<blockquote>
<b>Algorithm 3: Recursion</b>
<ul>
<li>Input: $d,r^2 \in \mathbb{N}$.</li>
<li>Output: $x \in B_\mathbb{Z}(d,r)$ uniformly random. </li>
</ul>
<ol>
<li>Allocate a lookup table for $f : \{0, \cdots, d\} \times \{0, \cdots, r^2\} \to \mathbb{N}$.</li>
<li>Initialize $f(0,s)=1$ for $s = 0 \cdots r^2$.</li>
<li>For $m = 1 \cdots d$:
<li>&nbsp;&nbsp;&nbsp;For $s = 0 \cdots r^2$:</li>
<li>&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;&nbsp;Set $f(m,s) = \sum_{a = \lceil - \sqrt{s} \rceil}^{\lfloor \sqrt{s} \rfloor} f(m-1, s-a^2 ).$</li>
<li>Set $s_1=r^2$ and allocate space for $x \in \{\lceil -r \rceil, \cdots ,\lfloor r \rfloor\}^d$.</li>
<li>For $i = 1 \cdots d$:</li>
<li>&nbsp;&nbsp;&nbsp;Sample $U_i \in \{1, \cdots, ~ f(d+1-i,s_i)\}$ uniformly at random.</li>
<li>&nbsp;&nbsp;&nbsp;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.$$</li>
<li>&nbsp;&nbsp;&nbsp;Set $s_{i+1}=s_i-x_i^2$.</li>
<li>Return $x \in B_\mathbb{Z}(d,r)$.</li>
</ol>
</blockquote>

<p>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)$.</p>

<hr>

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

<hr>

<blockquote>
<b>Algorithm 1: Rejection Sampling</b>
<ol>
<li>Sample $U_1, \cdots, U_d \in [-1,1]$ independently and uniformly at random.</li>
<li>Let $x = (U_1, \cdots, U_d)$.</li>
<li>If $\|x\|_2 \leq 1$, return $x$; else GOTO 1.</li>
</ol>
</blockquote>

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

<hr>

<blockquote>
<b>Algorithm 2: Scaled Spherical Sample</b>
<ol>
<li>Sample $G_1, \cdots, G_d \in \mathbb{R}$ independently from a standard Gaussian distribution ($\mathcal{N}(0,1)$).</li>
<li>Let $x = (G_1, \cdots, G_d)$ and $y=x/\|x\|_2$.</li>
<li>Sample $U \in [0,1]$ uniformly at random.</li>
<li>Return $z=y \cdot U^{1/d}$.</li>
</ol>
</blockquote> 

<p>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$}}.$$
<a href="https://en.wikipedia.org/wiki/Volume_of_an_n-ball" rel="nofollow noreferrer">We have</a> $$\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 <a href="https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform" rel="nofollow noreferrer">Box-Muller transform</a>.</p>

CC BY-SA 3.0
2016-02-07T19:48:16.920Z
I like Peter and your suggestions. There may be a little bias that I will check it
Bin Fu
CC BY-SA 3.0
2016-02-07T19:52:27.437Z
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.
Thomas Steinke
CC BY-SA 3.0
2016-02-07T21:20:00.447Z
It seems that function f(.) depends on the center if the original ball center
Bin Fu
CC BY-SA 3.0
2016-02-07T21:21:01.753Z
is not at a lattice point. Otherwise, it is uniform.
Bin Fu
CC BY-SA 3.0
2016-02-07T21:23:10.160Z
I am assuming the ball is centered at 0. Is this not the case?
Thomas Steinke
CC BY-SA 3.0
2016-02-07T21:25:17.620Z
@BinFu Can you please update the question to precisely define what distribution you want to sample from?
Thomas Steinke
CC BY-SA 3.0
2016-02-07T21:45:00.170Z
I look for the uniform random sampling for lattice points inside a d-dimension ball. If the center is at 0, the algorithm that you described looks uniform. There is a little bias when the center is not at a lattice point.  Overall, I am very happy for what I learnt here. Thanks
Bin Fu
CC BY-SA 3.0
2016-02-07T21:54:02.297Z
Algorithm 3 should generalize to sampling from a ball with a different center.
Thomas Steinke
CC BY-SA 3.0
2016-02-07T15:29:03.240Z
2016-02-07T20:22:44.677Z
2017-04-13T12:32:28.167Z
Thomas Steinke
<p>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.</p>

CC BY-SA 3.0
2016-02-07T06:55:29.197Z
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
Bin Fu
CC BY-SA 3.0
2016-02-07T06:56:45.057Z
Want to see a polynomial time solution.
Bin Fu
CC BY-SA 3.0
2016-02-07T06:49:30.587Z
2016-02-07T06:49:30.587Z
Jaeba
<p>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.</p>

<p>Thanks</p>

<p>Bin Fu</p>

CC BY-SA 3.0
2016-02-07T06:07:34.720Z
What do you mean by "grid point"? Do you just want to sample uniformly from the L2 ball?
Thomas Steinke
CC BY-SA 3.0
2016-02-07T08:20:34.417Z
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.
Rahul Savani
CC BY-SA 3.0
2016-02-07T14:21:00.237Z
What distribution? Uniform?
usul
CC BY-SA 3.0
2016-02-07T14:51:51.437Z
Yes, I am interested in lattice point as you pointed out. Your algorithm makes some sense.  We need to control time and bias.
Bin Fu
CC BY-SA 3.0
2016-02-07T15:04:47.960Z
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.
Bin Fu
CC BY-SA 3.0
2016-02-07T15:24:11.160Z
@BinFu: the easiest way to improve may be to not use a cube, but to use a slightly larger ball. If the radius is too small, it's still exponential. See my answer.
Peter Shor 
CC BY-SA 3.0
2016-02-07T21:40:05.293Z
@BinFu Please update your question to be precise about what you want. You have not clearly defined the distribution you want to sample from.
Thomas Steinke
<randomized-algorithms>
Random grid point in a d-dimensional ball
Bin Fu
# Random grid point in a d-dimensional ball

(Asked by: username_0 on 2016-02-07)

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

**Tags:** randomized-algorithms
> **Comments:**
> **username_1**: What do you mean by "grid point"? Do you just want to sample uniformly from the L2 ball?  
> **username_2**: 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.  
> **username_3**: What distribution? Uniform?  
> **username_0**: Yes, I am interested in lattice point as you pointed out. Your algorithm makes some sense.  We need to control time and bias.  
> **username_0**: 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.  
> **username_4**: @BinFu: the easiest way to improve may be to not use a cube, but to use a slightly larger ball. If the radius is too small, it's still exponential. See my answer.  
> **username_1**: @BinFu Please update your question to be precise about what you want. You have not clearly defined the distribution you want to sample from.  

---

## 3 answers

### Answer by username_4 on 2016-02-07. Score: 3

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 by username_1 on 2016-02-07. Score: 1

**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](https://cstheory.stackexchange.com/a/33755/2022) 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$:- For $s = 0 \cdots r^2$:
>    - Set $f(m,s) = \sum\_{a = \lceil - \sqrt{s} \rceil}^{\lfloor \sqrt{s} \rfloor} f(m-1, s-a^2 ).$
>    - Set $s\_1=r^2$ and allocate space for $x \in \{\lceil -r \rceil, \cdots ,\lfloor r \rfloor\}^d$.
>    - For $i = 1 \cdots d$:
>    - Sample $U\_i \in \{1, \cdots, ~ f(d+1-i,s\_i)\}$ uniformly at random.
>    - 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.$$
>    - Set $s\_{i+1}=s\_i-x\_i^2$.
>    - 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](https://en.wikipedia.org/wiki/Volume_of_an_n-ball) $$\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](https://en.wikipedia.org/wiki/Box%E2%80%93Muller_transform).

> **Comments:**
> **username_0**: I like Peter and your suggestions. There may be a little bias that I will check it  
> **username_1**: 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.  
> **username_0**: It seems that function f(.) depends on the center if the original ball center  
> **username_0**: is not at a lattice point. Otherwise, it is uniform.  
> **username_1**: I am assuming the ball is centered at 0. Is this not the case?  
> **username_1**: @BinFu Can you please update the question to precisely define what distribution you want to sample from?  
> **username_0**: I look for the uniform random sampling for lattice points inside a d-dimension ball. If the center is at 0, the algorithm that you described looks uniform. There is a little bias when the center is not at a lattice point.  Overall, I am very happy for what I learnt here. Thanks  
> **username_1**: Algorithm 3 should generalize to sampling from a ball with a different center.  

---

### Answer by username_5 on 2016-02-07. Score: 0

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.

> **Comments:**
> **username_0**: 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  
> **username_0**: Want to see a polynomial time solution.  

---

