Editing Low-discrepancy sequence (section)

==Construction of low-discrepancy sequences==

Because any distribution of random numbers can be mapped onto a uniform distribution, and quasirandom numbers are mapped in the same way, this article only concerns generation of quasirandom numbers on a multidimensional uniform distribution.

There are constructions of sequences known such that
:<math>
D_N^{*}(x_1,\ldots,x_N)\leq C\frac{(\ln N)^{s}}{N}.
</math>
where <math>C</math> is a certain constant, depending on the sequence. After Conjecture 2, these sequences are believed to have the best possible order of convergence. Examples below are the [[van der Corput sequence]], the [[Halton sequence]]s, and the [[Sobol sequence|Sobol’ sequence]]s. One general limitation is that construction methods can usually only guarantee the order of convergence. Practically, low discrepancy can be only achieved if <math>N</math> is large enough, and for large given s this minimum <math>N</math> can be very large. This means running a Monte-Carlo analysis with e.g. <math>s=20</math> variables and <math>N=1000</math> points from a low-discrepancy sequence generator may offer only a very minor accuracy improvement {{Citation needed|date=July 2023}}.

===Random numbers===

Sequences of quasirandom numbers can be generated from random numbers by imposing a negative correlation on those random numbers. One way to do this is to start with a set of random numbers <math>r_i</math> on <math>[0,0.5)</math> and construct quasirandom numbers <math>s_i</math> which are uniform on <math>[0,1)</math> using:

<math>s_i = r_i</math> for <math>i</math> odd and <math>s_i = 0.5 + r_i</math> for <math>i</math> even.

A second way to do it with the starting random numbers is to construct a random walk with offset 0.5 as in:

: <math>s_i = s_{i-1} + 0.5+ r_i \pmod 1. \, </math>

That is, take the previous quasirandom number, add 0.5 and the random number, and take the result [[modular arithmetic|modulo]]&nbsp;1.

For more than one dimension, [[Latin squares]] of the appropriate dimension can be used to provide offsets to ensure that the whole domain is covered evenly.

[[Image:Subrandom 2D.gif|thumb|270px|right|Coverage of the unit square. Left for additive quasirandom numbers with ''c''&nbsp;=&nbsp;0.5545497...,&nbsp;0.308517... Right for random numbers. From top to bottom. 10, 100, 1000, 10000 points.]]

===Additive recurrence===

For any irrational <math>\alpha</math>, the sequence

: <math>s_n = \{s_0 + n\alpha\}</math>

has discrepancy tending to <math>1/N</math>. Note that the sequence can be defined recursively by
: <math>s_{n+1} = (s_n + \alpha)\bmod 1 \;.</math>

A good value of <math>\alpha</math> gives lower discrepancy than a sequence of independent uniform random numbers.

The discrepancy can be bounded by the [[approximation exponent]] of <math>\alpha</math>. If the approximation exponent is <math>\mu</math>, then for any <math>\varepsilon>0</math>, the following bound holds:<ref name="kn05">{{Harvnb|Kuipers|Niederreiter|2005|p=123}}</ref>

: <math> D_N((s_n)) = O_\varepsilon (N^{-1/(\mu-1)+\varepsilon}).</math>

By the [[Thue–Siegel–Roth theorem]], the approximation exponent of any irrational algebraic number is 2, giving a bound of <math>N^{-1+\varepsilon}</math> above.

The recurrence relation above is similar to the recurrence relation used by a [[linear congruential generator]], a poor-quality pseudorandom number generator:<ref>{{Cite book|first=Donald E. |last= Knuth |title=[[The Art of Computer Programming]] |volume=2 |chapter=Chapter 3 – Random Numbers}}</ref>

: <math>r_i = (a r_{i-1} + c) \bmod m</math>

For the low discrepancy additive recurrence above, ''a'' and ''m'' are chosen to be 1.  Note, however, that this will not generate independent random numbers, so should not be used for purposes requiring independence.  

The value of <math>c</math> with lowest discrepancy is the fractional part of the [[golden ratio]]:<ref>
{{Cite web|first=Malte |last=Skarupke|url=https://probablydance.com/2018/06/16/fibonacci-hashing-the-optimization-that-the-world-forgot-or-a-better-alternative-to-integer-modulo/ |title=Fibonacci Hashing: The Optimization that the World Forgot|date=16 June 2018 |quote=One property of the Golden Ratio is that you can use it to subdivide any range roughly evenly ... if you don’t know ahead of time how many steps you’re going to take}}
</ref>

: <math>c = \frac{\sqrt{5}-1}{2} = \varphi - 1 \approx 0.618034.</math>

Another value that is nearly as good is the fractional part of the [[silver ratio]], which is the fractional part of the square root of 2:

: <math>c = \sqrt{2}-1 \approx 0.414214. \, </math>

In more than one dimension, separate quasirandom numbers are needed for each dimension. A convenient set of values that are used, is the square roots of [[primes]] from two up, all taken modulo 1:

: <math>c = \sqrt{2}, \sqrt{3}, \sqrt{5}, \sqrt{7}, \sqrt{11}, \ldots \, </math>

However, a set of values based on the generalised golden ratio has been shown to produce more evenly distributed points. <ref>{{Cite web|first=Martin |last=Roberts|year= 2018 |url=https://web.archive.org/web/20250301162105/https://extremelearning.com.au/unreasonable-effectiveness-of-quasirandom-sequences/ |title=The Unreasonable Effectiveness of Quasirandom Sequences}}
</ref>

The [[list of pseudorandom number generators]] lists methods for generating independent pseudorandom numbers.
''Note'': In few dimensions, recursive recurrence leads to uniform sets of good quality, but for larger <math>s</math> (like <math> s>8 </math>) other point set generators can offer much lower discrepancies.

===van der Corput sequence===
{{Main article|van der Corput sequence}}

Let

:<math>
n=\sum_{k=0}^{L-1}d_k(n)b^k
</math>

be the <math>b</math>-ary representation of the positive integer <math>n \geq 1</math>, i.e. <math>0 \leq d_k(n) < b</math>. Set

:<math>
g_b(n)=\sum_{k=0}^{L-1}d_k(n)b^{-k-1}.
</math>

Then there is a constant <math>C</math> depending only on <math>b</math> such that <math>(g_b(n))_{n \geq 1}</math> satisfies

:<math>
D^*_N(g_b(1),\dots,g_b(N))\leq C\frac{\log N}{N},
</math>

where <math>D^*_N</math> is the  
'''[[#Definition of discrepancy|star discrepancy]]'''.

===Halton sequence===
[[Image:Halton sequence 2D.svg|thumb|right|First 256 points of the (2,3) Halton sequence]]
{{Main article|Halton sequence}}

The Halton sequence is a natural generalization of the van der Corput sequence to higher dimensions. Let ''s'' be an arbitrary dimension and ''b''<sub>1</sub>, ..., ''b''<sub>''s''</sub> be arbitrary [[coprime]] integers greater than 1. Define

:<math>
x(n)=(g_{b_1}(n),\dots,g_{b_s}(n)).
</math>

Then there is a constant ''C'' depending only on ''b''<sub>1</sub>, ..., ''b''<sub>''s''</sub>, such that sequence {''x''(''n'')}<sub>''n''≥1</sub> is a ''s''-dimensional sequence with

:<math>
D^*_N(x(1),\dots,x(N))\leq C'\frac{(\log N)^s}{N}.
</math>

===Hammersley set===

[[File:Hammersley set 2D.svg|thumb|right|2D Hammersley set of size 256]]
Let <math>b_1,\ldots,b_{s-1}</math> be [[coprime]] positive integers greater than 1. For given <math>s</math> and <math>N</math>, the <math>s</math>-dimensional [[John Hammersley|Hammersley]] set of size <math>N</math> is defined by<ref name="HammersleyHandscomb1964">{{cite book|title=Monte Carlo Methods|last1=Hammersley|first1=J. M.|last2=Handscomb|first2=D. C.|year=1964|doi=10.1007/978-94-009-5819-7|isbn=978-94-009-5821-0 }}</ref>

:<math>
x(n)=\left(g_{b_1}(n),\dots,g_{b_{s-1}}(n),\frac{n}{N}\right)
</math>

for <math>n = 1, \ldots, N</math>. Then

:<math>
D^*_N(x(1),\dots,x(N))\leq C\frac{(\log N)^{s-1}}{N}
</math>

where <math>C</math> is a constant depending only on <math>b_1, \ldots, b_{s-1}</math>.

''Note'': The formulas show that the Hammersley set is actually the Halton sequence, but we get one more dimension for free by adding a linear sweep. This is only possible if  <math>N</math> is known upfront. A linear set is also the set with lowest possible one-dimensional discrepancy in general. Unfortunately, for higher dimensions, no such "discrepancy record sets" are known. For <math>s=2</math>, most low-discrepancy point set generators deliver at least near-optimum discrepancies.

===Sobol sequence===
{{Main article|Sobol sequence}}

The Antonov–Saleev variant of the Sobol’ sequence generates numbers between zero and one directly as binary fractions of length <math>w,</math> from a set of <math>w</math> special binary fractions, <math>V_i, i = 1, 2, \dots, w</math> called direction numbers. The bits of the [[Gray code]] of <math>i</math>, <math>G(i)</math>, are used to select direction numbers. To get the Sobol’ sequence value <math>s_i</math> take the [[exclusive or]] of the binary value of the Gray code of <math>i</math> with the appropriate direction number. The number of dimensions required affects the choice of <math>V_i</math>.

=== Poisson disk sampling ===
{{main article|Supersampling#Poisson disk}}
[[Supersampling#Poisson disk|Poisson disk sampling]] is popular in video games to rapidly place objects in a way that appears random-looking
but guarantees that every two points are separated by at least the specified minimum distance.<ref>
Herman Tulleken.
{{Cite magazine|url=http://devmag.org.za/2009/05/03/poisson-disk-sampling/ |title=Poisson Disk Sampling |magazine=Dev.Mag | issue= 21 |date= March 2008| first=Herman |last=Tulleken |pages=21–25 }}
</ref> This does not guarantee low discrepancy (as e. g. Sobol’), but at least a significantly lower discrepancy than pure random sampling. The goal of these sampling patterns is based on frequency analysis rather than discrepancy, a type of so-called "blue noise" patterns.