17 ms·
The Unreasonable Effectiveness of Quasirandom Sequences
- charlesism 8y agoNeat! Looks like a pretty decent way to generate colored noise. Wish there were an audio demo.
- extremelearning 8y agoOnly blue noise... I have absolutely no idea what it would sound like. So that means I’ll definitely try it later tonight and update the blog post. Tx
- charlesism 8y agoJudging by the nicely spaced points in the first animation on your web page, probably pretty good.
- extremelearning 8y agoI have now inserted some sound demonstrations into my blog post! You can now listen to what 1-dimensional (mon) and 2-dimensional (stereo) quasirandom sound might be like. I will leave it to you to decide which ones you prefer, and maybe even think why some of the quasirandom sequences are more pleasant than others. For those who have already visited the site, you may need to clear your browser cache to see/hear this update.
- charlesism 8y agoThanks! Interesting. At some point I want to play around with this myself. I'm curious to see what it would sound like to combine the x and y values together (eg: if the frequency is 10 samples between each x, randomly choose between x=0 to x=10 for the first, between x=10 to x=20 for the second). That would probably make it sound less like an FM tone, and more like hissing. I'm not much for math, but a lot of tutorials on colored noise are based on multiple filters. If this algorithm works, it would just need a single low-pass filter for smoothing.
- Tomte 8y ago[I failed to reply to the author, but I guess you'll see it here, as well] Please remove the plea for upvotes from your site. The mods usually react harshly to this and bury the submission. Your article doesn't need that to stay some time on the front page, anyway. It's high-quality and people have recognized it as such.
- extremelearning 8y agoThanks for this advice. I am new to blogging and certainly new to HN front page etiquette. I think my excitement got the better of me! I have now taken it off.
- Tomte 8y agoGreat, thanks!
- PavlikPaja 8y agoPlease note the article only shows correctly in Chrome.
- PavlikPaja 8y agoIt seems to be fixed now.
- deleted 8y ago[deleted]
- strainer 8y agoI would like to be able to make this sequence, it looks nice. ~~Unfortunately the included python code is not valid~~. ( * ) It will take me quite a while to grok the math discussion well enough to implement it. * - pasted it wrong.
- extremelearning 8y agoThat's strange. To double-check that i typed it into my blog correctly, I copy-and-pasted the code from the post into my python IDE and it compiles and runs properly in my environment. For what it is worth, I am using Python 2.7 on a 64-bit windows machine. Maybe someone else reading has some ideas? You should notice that after my example code, I have included two other people's code demos so that you can see the same algorithm from a different perspective. Maybe these will help you understand what's required. In summary, for 2 dimensions, x and y coordinates of the n-th term (n = 1,2,3,....) are defined as: g = 1.32471795724474602596090885447809 a1 = 1.0/g a2 = 1.0/(g*g) x[n] = (0.5+a1*n) %1 y[n] = (0.5+a2*n) %1 where %1 is the mod 1 operator and takes the fractional part of the argument. Hope that helps.
- strainer 8y agoI pasted it wrong blush And, Thankyou that seems as clear as can be, when I get my head screwed on straight I will have a better read through your fine insights of the workings. I will be trying include your sequence in my quirky collection of random helpers for Javascript, illustrated here. http://strainer.github.io/Fdrandom.js/ http://strainer.github.io/Fdrandom.js/ And will be sure to attribute you and message when I get it in.
- extremelearning 8y agoThat's great. I'm glad you got it working. ;)
- bzzzzzt 8y agoThis looks similar to the 2D "random number generators" at http://pippin.gimp.org/a_dither/ http://pippin.gimp.org/a_dither/ but is even more evenly distributed.
- extremelearning 8y agoYes. I think that using the dither matrix based on the R_2 sequence should be competitive with any of the the dither masks mentioned on that link. Thus in situations like video, and/or ray tracing where computing speed is of the essence, maybe a better mask such as the R2 dither mask would offer an improvement on existing methods. However, in terms of end-quality I am not aware of any dither masks that are anywhere near as good as the state-of-the-art error diffusion algorithms. Thus, for a small number of images such as our favourite animated gifs, I would still be recommending an error-diffusion method. :)
- targafarian 8y agoI'm curious why Sobol outperforms the R-sequence for numerical integration of an 8-dimensional Gaussian and whether that result holds in lower dimensions (and whether that's tied to the infintely-long-tailed nature of the Gaussian; maybe there's a bias in Sobol that is beneficial in this particular case...?).
- extremelearning 8y agoSobol's sequence is really quite amazing. However, there are three reasons why it doesn't get as much attention as it probably deserves. First is that the maths is far more complex than that of other sequences. Here is a the easiest description I have seen so far [1]. Secondly this complexity as well as the non-trivial requirement to carefully select the basis parameters means that the computing side is very very complex. To get an indication of its complexity, the code for Sobol generation in Mathematica and WolframAlpha is not done by Wolfram itself but rather "comes courtesy of the Intel MKL libraries,... Specifics of the implementation, such as choice of initial values, are not documented. Evaluation of this implementation in terms of discrepancy, projections, and performance in application remains for future work." [also 1] Thirdly, as defined by d* discrepancy, ( which is by far the most common formal method of quatifying discrepancy) the Sobol sequence does not have an asymptotic discrepancy as low as many of other sequences. Despite this last point, it has been found that in practice Sobol generally yields as good or better performance when used for numerical integrations. I personally have a suspicion that this is because of the special property, called "Property A" [2] that all Sobol point sequences possess. However, I will let people much smarter than me comment on whether this belief is justified or not... In one of my other posts [3], Sobol sequences outperform all of the other low discrepancy sequence (including my new R2 sequence) hands down. Finally, please note that for brevity and clarity, my blog only shows one example for most topics. I will leave it to other authors, or future posts to more rigourously show how consistent the comparisons are.... [1] http://www.mathematica-journal.com/2011/12/a-toolbox-for-quasirandom-simulation/ http://www.mathematica-journal.com/2011/12/a-toolbox-for-qua... [2] http://www.andreasaltelli.eu/file/repository/HD_SobolGenerator.pdf http://www.andreasaltelli.eu/file/repository/HD_SobolGenerat... [3] http://extremelearning.com.au/how-to-generate-a-sequence-of-evenly-balanced-combinations/ http://extremelearning.com.au/how-to-generate-a-sequence-of-...
- deleted 8y ago[deleted]
- tlb 8y agoI'm going to try this instead of Halton sequences in my current project. One detail of the implementation that may be worth improving: z[i] = (seed + alpha*(i+1)) %1 When i gets large, say 2^30, you're losing 30 bits of precision in the %1, so (with IEEE doubles) it can only take one of 2^23 values. The usual implementation of Halton sequences get slower with large values of i, but doesn't suffer from quantization.
- extremelearning 8y agoThanks tlb. You're totally right. I wrote the example code to match the notation used in the post which is a reflection of my maths background. I suspect that from a programming perspective using the recurrence relation: z[i+1] = (z[i] + alpha) %1 is better practice in terms of both speed and accuracy. I have updated my post to make this clearer. Thanks.
- extremelearning 8y agoAuthor here. OMG! Someone submitted my post here and it's got to the front page!! Happy to try to answer any questions that people have. Enjoy!
- leni536 8y agoI also played with improving the Halton sequence. There is an other generalization of the Golden ratio, namely the metallic ratios. Unfortunately using them performs worse than the Halton sequence.
- nestorD 8y agoVery nice blog (I subscribed to the newsletter). I am looking forward to the "quasirandom stochastic descent" and "extreme learning machines" sections.
- SeanLuke 8y agoForgive me if you make this clear in the text (I may have missed it), but you have various pull quotes which say: "The new R2 sequence is the only 2-dimensional low discrepancy quasirandom sequence where..." Are you saying that it's the only known sequence? Or the only one possible? If the second, do you have formal proofs for those claims?
- extremelearning 8y agoI somewhat stumbled on this during my efforts to improve some machine learning architecture I was working on, so unfortunately no formal proofs as yet. This is one reason why I went down the route of technical blog post rather than more academic routes which would no doubt require more rigour and precise language. Although in many cases I have reason to strongly believe that it is not possible for other sequences to have similar properties, until such proofs are established, I agree that in a technical paper it would be prudent for me to preface most of the quotes with “only known...”
- jacobolus 8y agoOne thing I wondered is if there’s a way to tweak this (e.g. by skipping values in one direction) to have better region filling (with respect to Euclidean distance) when the aspect ratio is not square. Try sliding the top slider at https://beta.observablehq.com/@jrus/plastic-sequence https://beta.observablehq.com/@jrus/plastic-sequence to see what I mean.
- extremelearning 8y agoUnfortunately, I haven’t found a nice solution to this problem yet. For those requiring quasirandom Monte Carlo integration over a say a 2x1 rectangular grid there is always the option of scaling the unit square points by 2 and then rejecting half of the points. Not super elegant and not very efficient in high dimensions, but at least it is unbiased and maintains the low discrepancy characteristics.
- 8y ago
- mturmon 8y agoAnyone who does Monte Carlo simulations would at least be interested in the notion of Quasi Monte Carlo, in which points are sampled "more evenly than random." As noted in OP, this can lead to O(1/n) convergence of averages to underlying target values, rather than the expected O(1/sqrt(n)) convergence ("n" is the number of probes of the underlying system). This is shown, strikingly, in Fig. 2 of the OP. There is a conference, tutorial papers etc. E.g.: http://mcqmc2018.inria.fr http://mcqmc2018.inria.fr
- mlthoughts2018 8y agoI work professionally in large-scale mcmc projects and don’t find quasi mc to be much more than a curiosity in my work, because integration (as opposed to drawing a rich set of posterior samples) is virtually never the goal. To do anything of use, you have to draw samples and then simulate complex outcomes, very few of which reduce to simple averages over the drawn samples or functions of the drawn samples. And since you need to draw samples, you might as well use those samples for any parts that do happen to require integration, rather than a quasi mc calculation, except in extremely toy-problem situations where the convergence rate matters disproportionately for that smallset of outcomes. I agree that for cases when you just want to evaluate an integral it could be useful. I have never encountered a use case when anyone just wanted to calculate an integral, as opposed to also generating posterior uncertainty metrics, posterior test statistics for posterior predictive checking, posterior diagnostics like ordinal statistics among the posterior samples. I’m sure outside of stats, there must be use cases. Just adding a counterpoint to the idea that quasi mc should always be interesting to practitioners. For a lot of people who work in mcmc methods, quasi mc is just not interesting and generally speaking could never be.
- mturmon 8y agoI agree with all of this. QMC is more relevant to Monte Carlo, but MCMC is its own thing as far as I know. There is some work on combining the two, but it doesn’t seem mature. If you don’t mind: what’s the general application domain for your MCMC‘s? For me, it’s inference of (usually) spatial fields for various science data sets, like Earth or solar imagery or atmospheric composition.
- infogulch 8y ago> The methods of creating fully deterministic low discrepancy quasirandom sequences in one dimension are extremely well studied, and generally solved. Are you referring to something like quadratic residue [1]? It's basically that the sequence x^2 mod P where P is prime appears to be somewhat random, and it goes through the whole sequence before repeating. But it's obviously finite, and it's not perfect [2]. I'm mostly just curious how it relates. [1]: https://en.wikipedia.org/wiki/Quadratic_residue https://en.wikipedia.org/wiki/Quadratic_residue [2]: https://arxiv.org/ftp/arxiv/papers/1612/1612.05852.pdf https://arxiv.org/ftp/arxiv/papers/1612/1612.05852.pdf
- extremelearning 8y agoMy intent when I wrote that phrase was two-fold. Firstly to emphasise the large body of work that was done by Weyl, Kronecker, Hurwitz et al relating to the equidistribtuion theorem, badly approximately numbers and diophantine approximation. Armed with this knowledge it is easily proved that the recurrrence relation using the golden ratio mod 1 is optimal using almost any measure. Although some other readers on HN may be able to see a connection between one dimensional quasirandom sequences and quadratic residue, unfortunately I can’t immediately see an explicit one. Sorry! My second intent was to contrast this with how little is provably known for higher dimensions. To my knowledge, thanks to the phenomenal work by Halton, Sobol, Niederreiter et al, we know that many of the contemporary sequences are optimal in the limiting Big O sense, but there are no proofs even for d=2 for optimality for any of the sequences for general finite values of n.
- leni536 8y agoVery cool! My only gripe of low discrepancy sequences is how they are used in Quasi Monte-Carlo integration: 1/n \sum_1^N f(a_n), so each point gets equal weight. One can generalize this that every point gets it's own weight: \sum_1^N w_{N,n} f(a_n). I have never seen this generalization elsewhere. I only study Quasi Monte-Carlo as a hobby, maybe I just missed it. Once I had significantly improved discrepancies (close to theoretical minimum) with a modified Van der Corput sequence and linearly diminishing weights. It's all in one dimension though. My other idea for two dimensions were to use the Hilbert-curve (the infinite limit) to map a one-dimensional sequence to a square area, but AFAIK this was done by others before.
- extremelearning 8y agoYes, using an open (infinite) one dimensional low-discrepancy sequence in conjunction with a space-filling curve such as the Hilbert-curve has been explored. For example by Owen [1] Niederreiter himself has a beautiful paper that includes a section on weighted low discrepancy sequence [2] where he also cites several references within it that might be of interest to you. [1] http://statweb.stanford.edu/~owen/reports/extgrid.pdf http://statweb.stanford.edu/~owen/reports/extgrid.pdf [2] http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.614.2202&rep=rep1&type=pdf http://citeseerx.ist.psu.edu/viewdoc/download?doi=10.1.1.614...