You signed in with another tab or window. Reload to refresh your session.You signed out in another tab or window. Reload to refresh your session.You switched accounts on another tab or window. Reload to refresh your session.Dismiss alert
Copy file name to clipboardexpand all lines: contents/box_muller/box_muller.md
+10-173
Original file line number
Diff line number
Diff line change
@@ -1,6 +1,6 @@
1
1
# The Box—Muller Transform
2
2
3
-
The Box—Muller transform holds a special pace in my heart as it was the first method I ever had to implement for my own research.
3
+
The Box—Muller transform holds a special place in my heart as it was the first method I ever had to implement for my own research.
4
4
5
5
The purpose of this transformation is simple.
6
6
It takes a uniform (probably random) distribution and turns it into a Gaussian one.
@@ -11,8 +11,8 @@ It takes a uniform (probably random) distribution and turns it into a Gaussian o
11
11
12
12
That's it.
13
13
14
-
It was originally developed by George Box (yes, Box is his last name) and Mervin Muller in 1958 and is one of the most common methods to create a random, Gaussian distributions of points {{ "box1958" | cite }}.
15
-
It's particularly useful when initializing a distribution of particles for a physical, N-body simulation.
14
+
It was originally developed by George Box (yes, Box is his last name) and Mervin Muller in 1958 and is one of the most common methods to create a random, Gaussian distribution of points {{ "box1958" | cite }}.
15
+
It's particularly useful when initializing a set of particles for a physical, N-body simulation.
16
16
This chapter will be divided into a few subsections:
17
17
18
18
1. How to initialize the Box—Muller transform
@@ -46,7 +46,7 @@ Note that there is an inherent limitation with this method in that it only works
46
46
Because of this, we decided to round $$n$$ up to the nearest square to make a nice grid.
47
47
It's not the cleanest implementation, but the grid will mainly be used for debugging anyway, so it's OK to be a *little* messy here.
48
48
49
-
The real star of the show here is the uniform random distribution, which can be generated like this:
49
+
The real star of the show is the uniform random distribution, which can be generated like this:
50
50
51
51
{% method %}
52
52
{% sample lang="jl" %}
@@ -65,19 +65,18 @@ Good question!
65
65
The easiest way is to plot a histogram of a super large number of points.
66
66
If the random distribution is uniform, then all the bins should be roughly the same value.
67
67
The more points we have, the smaller the percent difference between the bins will be.
68
-
Here is a set of images for $$n=100$$, $$1,000$$, and $$10,000$$ all in one dimension:
68
+
Here is a set of images generated by `rand()`for $$n=100$$, $$1,000$$, and $$10,000$$ all in one dimension:
It is clear that the 10,000 case looks the most uniform.
75
+
It is clear that the $$10,000$$ case looks the most uniform.
76
76
Note that for two dimensions, the same logic applies, but we need to create separate histograms for the $$x$$ and $$y$$ coordinates.
77
77
78
78
Once this test is complete, we can be fairly sure that the function we are using to generate a random distribution is uniform and ready for the next step of the process: actually using the Box—Muller transform!
79
79
80
-
81
80
## How to use the Box—Muller transform in Cartesian coordinates
82
81
83
82
The (two dimensional) Cartesian version of the Box—Muller transform starts with two random input values ($$u_1$$ and $$u_2$$), both of which come from their own uniform distributions that are between $$0$$ and $$1$$.
@@ -106,10 +105,6 @@ y &= r\sin(\theta).
106
105
\end{align}
107
106
$$
108
107
109
-
This can be more easily visualized in the following animation:
110
-
111
-
ADD ANIMATION
112
-
113
108
Finally, in order to specify the size and shape of the generated Gaussian distribution, we can use the standard deviation, $$\sigma$$, and the mean, $$\mu$$, like so:
114
109
115
110
$$
@@ -123,7 +118,7 @@ In general, this can all be written in code like so:
@@ -136,63 +131,11 @@ Which produces the following output
136
131
Note that we have written the code to work on a single set of input values, but it could also be written to read in the entire distribution of points all at once.
137
132
As this particular technique is usually implemented in parallel, it's up to you to decided which is the fastest for your own individual use-case.
138
133
139
-
For now, I would like to briefly touch on how to abstract this to different dimensions.
140
-
141
-
### What about 1D?
142
-
143
-
If we would like to use the Box—Muller transform in one dimension, things are relatively easy, just ignore $$u_2$$ and $$\theta$$, such that
144
-
145
-
$$
146
-
z = \sigma r + \mu = \sigma\sqrt{-2\ln(u)} + \mu.
147
-
$$
148
-
149
-
There shouldn't be any funny business here.
150
-
We will leave the code modifications to the reader, but show the results:
151
-
152
-
ADD IMAGE
153
-
154
-
### OK, what about for 3D?
155
-
156
-
For three dimensions, we need to use spherical coordinates and one additional uniform distribution ($$u_3$$) to create a third point $$z_3$$.
157
-
Also, the additional angular value of $$\phi$$ only stretches to $$\pi$$, so the appropriate transforms would be
158
-
159
-
$$
160
-
\begin{align}
161
-
r &= \sqrt{-2\ln(u_1)} \\
162
-
\theta &= 2\pi u_2 \\
163
-
\phi &= \pi u_3.
164
-
\end{align}
165
-
$$
166
-
167
-
Then $$x,y,z$$ are
168
-
169
-
$$
170
-
\begin{align}
171
-
x &= r\cos(\theta)\sin(\phi) \\
172
-
y &= r\sin(\theta)\sin(\phi) \\
173
-
z &= r\cos(\phi),
174
-
\end{align}
175
-
$$
176
-
177
-
and
178
-
179
-
$$
180
-
\begin{align}
181
-
z_1 &= \sigma x + \mu \\
182
-
z_2 &= \sigma y + \mu \\
183
-
z_3 &= \sigma z + \mu.
184
-
\end{align}
185
-
$$
186
-
187
-
Similar to the previous subsection, we will leave the code modifications to the reader and show the results:
188
-
189
-
ADD IMAGE
190
-
191
134
At this stage, we have a good idea of how the transform works, but some people shy away from the Cartesian method in practice and instead opt for the polar form, which will discussed next!
192
135
193
136
## How to use the Box—Muller transform in polar coordinates
194
137
195
-
The Cartesian form of the Box—Muller transform is relatively intuitive
138
+
The Cartesian form of the Box—Muller transform is relatively intuitive.
196
139
The polar method is essentially the same, but without the costly $$\sin$$ and $$\cos$$ operations.
197
140
In this case, we use the input values to create an initial radial point (to be scaled later):
198
141
@@ -231,10 +174,6 @@ z_2 &= \sigma y + \mu.
231
174
\end{align}
232
175
$$
233
176
234
-
This can be more easily visualized in the following animation:
235
-
236
-
ADD ANIMATION
237
-
238
177
In code, this might look like this:
239
178
240
179
{% method %}
@@ -250,112 +189,10 @@ This will produce the following output:
250
189
</p>
251
190
252
191
Again, this is ultimately the same as the Cartesian method, except that it:
253
-
1. Rejects points outside the unit circle (also called rejection sampling)
192
+
1. Rejects points in the initial distribution that are outside of the unit circle (also called rejection sampling)
254
193
2. Avoids costly $$\sin$$ and $$\cos$$ operations
255
194
256
-
Point 2 means that the polar method *should be* way faster than the Cartesian one, but rejection sampling is somewhat interesting in it's own right...
257
-
258
-
### Just how costly is rejection sampling anyway?
259
-
260
-
Let's imagine we want to have a final Gaussian with $$n$$ particles in it.
261
-
With the Cartesian Box—Muller method, this is easy: start the initial distribution(s) with $$n$$ particles and then do the transform.
262
-
Things *can* be just as easy with the Polar Box—Muller method as well, so long as we start with a uniformly distributed circle instead of a uniformly distributed square.
263
-
That is to say, so long as we do the rejection sampling before-hand, the Polar Box—Muller method will always be more efficient.
264
-
To be fair, there are methods to generate a uniform distribution of points within a circle without rejection sampling, but let's assume that we require rejection sampling for this example
265
-
266
-
This means that someone somehow needs to do the rejection sampling for the Polar method, which is sometimes a painful process.
267
-
This also means that the Box—Muller method can be used to teach some of the fundamentals of General-Purpose GPU computing.
268
-
Note that because of the specificity of this problem, all the code in this subsection will be in Julia and using the package KernelAbstractions.jl, which allows us to execute the same kernels on either CPU or GPU hardware depending on how we configure things.
269
-
270
-
Let's first consider the case where we do the rejection sampling as a part of the polar Box—Muller kernel instead of as a pre-processing step.
271
-
In this case, we can imagine 2 separate ways of writing our kernel:
272
-
1. With replacement: In this case, we *absolutely require* the final number of points in our Gaussian distribution to be $$n$$, so if we find a point outside of the unit circle while running the kernel, we will "re-roll" again for a new point that *is* within the circle.
273
-
2. Without replacement: This means that we will start with a uniform distribution of $$n$$ points, but end with a Gaussian of $$m < n$$ points. In this case, if we find a point outside of the unit circle while running the kernel, we just ignore it by setting the output values to NaNs (or something).
1. If we find a point outside of the unit circle, we have to continually look for new points until we *do* find one inside of the circle. This means that some threads might take literally forever to find a new point (if we are really unlucky).
282
-
2. To generate new points, we need to re-generate a uniform distribution, but what if our uniform distribution is not random? What if it's a grid (or something similar) instead? In this case, we really shouldn't look for a new point on the inside of the circle as all those points have already been accounted for.
283
-
3. The `rand()` function is kinda tricky on some parallel platforms (like GPUs) and might not work out of the box. In fact, the implementation shown above can only be run on the CPU.
284
-
285
-
OK, fine.
286
-
I don't think anyone expected a kernel with a `while` loop inside of it to be fast.
287
-
So what about a method without replacement?
288
-
Surely there is no problem if we just ignore the `while` loop altogether!
289
-
Well, the problem with this approach is a bit less straightforward, but first, code:
To start discussing why the kernel without replacement is *also* a bad idea, let's go back to the [Monte Carlo chapter](../monte_carlo/monte_carlo.md), where we calculated the value of $$\pi$$ by embedding it into a circle.
294
-
There, we found that the probability of a randomly chosen point falling within the unit circle to be $$\frac{\pi r^2}{(2r)^2} = \frac{pi}{4} \sim 78.54\%$$, shown in the visual below:
This means that a uniform distribution of points within a circle will reject $$\sim 21.46\%$$ of points on the square.
301
-
This also means that if we have a specific $$n$$ value we want for the final distribution, we will need $$\frac{1}{0.7853} \sim 1.273 \times$$ more input values on average!
302
-
303
-
No problem!
304
-
In this hypothetical case, we don't need *exactly*$$n$$ points, so we can just start the initial distributions with $$1.273 \times n$$ points, right?
305
-
306
-
Right.
307
-
That will work well on parallel CPU hardware, but on the GPU this will still have an issue.
308
-
309
-
On the GPU, computation is all done in parallel, but there is a minimum unit of parallelism called a *warp*.
310
-
The warp is the smallest number of threads that can execute something in parallel and is usually about 32.
311
-
This means that if an operation is queued, all 32 threads will do it at the same time.
312
-
If 16 threads need to execute something and the other 16 threads need to execute something else, this will lead to *warp divergence* where 2 actions need to be performed instead of 1:
In this image, every odd thread needs to perform the pink action, while the even threads need to perform the blue action.
319
-
This means that 2 separate parallel tasks will be performed, one for the even threads, another for the odd threads.
320
-
This means that if $$\ell$$ operations are queued, it could take $$\ell\times$$ as long for all the threads to do their work!
321
-
This is why `if` statements in a kernel can be dangerous!
322
-
If used improperly, they can cause certain threads in a warp to do different things!
323
-
324
-
A good question here is, "why can't we modify the warp so that it only uses odd threads or even threads?"
325
-
This would prevent the warp divergence entirely, but it comes at the cost of more complicated memory operations, which is an entirely different problem.
326
-
For the most part, even it our warp does diverge, it is a good idea to make sure that each thread is working on it's own array element without complex accessing patterns.
327
-
328
-
If we look at the above kernel, we are essentially asking $$78.53\%$$ of our threads to do something different than everyone else, and because we are usually inputting a uniform random distribution, this means that *most* warps will have to queue up 2 parallel actions instead of 1.
329
-
Now we need to pick our poison: slow $$\sin$$ and $$\cos$$ operations or warp divergence.
330
-
331
-
The only way to know which is better is to perform benchmarks, which we will show in a bit, but there is one final scenario we should consider: what about doing the rejection sampling as a pre-processing step?
332
-
This would mean that we pre-initialize the polar kernel with a uniform distribution of points in the unit circle.
333
-
This means no warp divergence, so we can get the best of both worlds, right?
334
-
335
-
Well, not exactly.
336
-
The polar Box—Muller method will definitely be faster, but again: someone somewhere needed to do rejection sampling and if we include that step into the process, things become complicated again.
337
-
The truth is that this pre-processing step is difficult to get right, so it might require a chapter in it's own right.
338
-
339
-
In many cases, it's worth spending a little time before-hand to make sure subsequent operations are fast, but in this case, we only have a single operation, not a set of operations.
340
-
The Box—Muller method will usually only be used once at the start of the simulation, which means that the pre-processing step of rejection sampling might end up being overkill.
341
-
342
-
No matter the case, benchmarks will show the true nature of what we are dealing with here:
| Polar with replacement |$$433.644 \pm 2.64$$ms | NA |
349
-
350
-
These were run with an Nvidia GTX 970 GPU and a Ryzen 3700X 16 core CPU.
351
-
Again, the code can be found [here](code/julia/performance.jl) and was written in Julia using the KernelAbstractions.jl package for parallelization.
352
-
For these benchmarks, we used Julia's inbuild benchmarking suite from BenchmarkTools, making sure to sync the GPU kernels with `CUDA.@sync`.
353
-
We also ran with $$4096^2$$ input points.
354
-
355
-
Here, we see an interesting divergence in the results.
356
-
On the CPU, the polar method is *always* faster, but on the GPU, both methods are comparable.
357
-
I believe this is the most important lesson to be learned from the Box—Muller method: sometimes, no matter how hard you try to optimize your code, different hardware can provide radically different results!
358
-
It's incredibly important to benchmark code to make sure it is actually as performant as you think it is!
195
+
Point 2 means that the polar method *should be* way faster than the Cartesian one, but rejection sampling is somewhat interesting in it's own right, which I have added in a [separate chapter](box_muller_rejection.md)
0 commit comments