Low-discrepancy sequence

Error in estimated kurtosis as a function of number of datapoints. 'Additive quasirandom' gives the maximum error when c = (√5 − 1)/2. 'Random' gives the average error over six runs of random numbers, where the average is taken to reduce the magnitude of the wild fluctuations

Quasirandom numbers have an advantage over pure random numbers in that they cover the domain of interest quickly and evenly.

Two useful applications are in finding the characteristic function of a probability density function, and in finding the derivative function of a deterministic function with a small amount of noise. Quasirandom numbers allow higher-order moments to be calculated to high accuracy very quickly.

Applications that don't involve sorting would be in finding the mean, standard deviation, skewness and kurtosis of a statistical distribution, and in finding the integral and global maxima and minima of difficult deterministic functions. Quasirandom numbers can also be used for providing starting points for deterministic algorithms that only work locally, such as Newton–Raphson iteration.

Quasirandom numbers can also be combined with search algorithms. A binary tree Quicksort-style algorithm ought to work exceptionally well because quasirandom numbers flatten the tree far better than random numbers, and the flatter the tree the faster the sorting. With a search algorithm, quasirandom numbers can be used to find the mode, median, confidence intervals and cumulative distribution of a statistical distribution, and all local minima and all solutions of deterministic functions.

The discrepancy of a set P = {x₁, ..., x_N} is defined, using Niederreiter's notation, as

${\displaystyle D_{N}(P)=\sup _{B\in J}\left|{\frac {A(B;P)}{N}}-\lambda _{s}(B)\right|}$

where λ_s is the s-dimensional Lebesgue measure, A(B;P) is the number of points in P that fall into B, and J is the set of s-dimensional intervals or boxes of the form

${\displaystyle \prod _{i=1}^{s}[a_{i},b_{i})=\{\mathbf {x} \in \mathbf {R} ^{s}:a_{i}\leq x_{i}<b_{i}\}\,}$

where ${\displaystyle 0\leq a_{i}<b_{i}\leq 1}$.

The star-discrepancy D^*_N(P) is defined similarly, except that the supremum is taken over the set J^* of rectangular boxes of the form

${\displaystyle \prod _{i=1}^{s}[0,u_{i})}$

where u_i is in the half-open interval [0, 1).

The two are related by

${\displaystyle D_{N}^{*}\leq D_{N}\leq 2^{s}D_{N}^{*}.\,}$

Note: With these definitions, discrepancy represents the worst-case or maximum point density deviation of a uniform set. However, also other error measures are meaningful, leading to other definitions and variation measures. For instance, L2 discrepancy or modified centered L2 discrepancy are also used intensively to compare the quality of uniform point sets. Both are much easier to calculate for large N and s.

Let Ī^s be the s-dimensional unit cube, Ī^s = [0, 1] × ... × [0, 1]. Let f have bounded variation V(f) on Ī^s in the sense of Hardy and Krause. Then for any x₁, ..., x_N in I^s = [0, 1) × ... × [0, 1),

${\displaystyle \left|{\frac {1}{N}}\sum _{i=1}^{N}f(x_{i})-\int _{{\bar {I}}^{s}}f(u)\,du\right|\leq V(f)\,D_{N}^{*}(x_{1},\ldots ,x_{N}).}$

The Koksma–Hlawka inequality is sharp in the following sense: For any point set {x₁,...,x_N} in I^s and any ${\displaystyle \varepsilon >0}$, there is a function f with bounded variation and V(f) = 1 such that

${\displaystyle \left|{\frac {1}{N}}\sum _{i=1}^{N}f(x_{i})-\int _{{\bar {I}}^{s}}f(u)\,du\right|>D_{N}^{*}(x_{1},\ldots ,x_{N})-\varepsilon .}$

Therefore, the quality of a numerical integration rule depends only on the discrepancy D^*_N(x₁,...,x_N).

Let ${\displaystyle D=\{1,2,\ldots ,d\}}$. For ${\displaystyle \emptyset \neq u\subseteq D}$ we write

${\displaystyle dx_{u}:=\prod _{j\in u}dx_{j}}$

and denote by ${\displaystyle (x_{u},1)}$ the point obtained from x by replacing the coordinates not in u by ${\displaystyle 1}$. Then

${\displaystyle {\frac {1}{N}}\sum _{i=1}^{N}f(x_{i})-\int _{{\bar {I}}^{s}}f(u)\,du=\sum _{\emptyset \neq u\subseteq D}(-1)^{|u|}\int _{[0,1]^{|u|}}\operatorname {disc} (x_{u},1){\frac {\partial ^{|u|}}{\partial x_{u}}}f(x_{u},1)\,dx_{u},}$

where ${\displaystyle \operatorname {disc} (z)={\frac {1}{N}}\sum _{i=1}^{N}\prod _{j=1}^{d}1_{[0,z_{j})}(x_{i,j})-\prod _{j=1}^{d}z_{i}}$ is the discrepancy function.

Applying the Cauchy–Schwarz inequality for integrals and sums to the Hlawka–Zaremba identity, we obtain an ${\displaystyle L^{2}}$ version of the Koksma–Hlawka inequality:

${\displaystyle \left|{\frac {1}{N}}\sum _{i=1}^{N}f(x_{i})-\int _{{\bar {I}}^{s}}f(u)\,du\right|\leq \|f\|_{d}\operatorname {disc} _{d}(\{t_{i}\}),}$

where

${\displaystyle \operatorname {disc} _{d}(\{t_{i}\})=\left(\sum _{\emptyset \neq u\subseteq D}\int _{[0,1]^{|u|}}\operatorname {disc} (x_{u},1)^{2}\,dx_{u}\right)^{1/2}}$

and

${\displaystyle \|f\|_{d}=\left(\sum _{u\subseteq D}\int _{[0,1]^{|u|}}\left|{\frac {\partial ^{|u|}}{\partial x_{u}}}f(x_{u},1)\right|^{2}dx_{u}\right)^{1/2}.}$

${\displaystyle L^{2}}$ discrepancy has a high practical importance because fast explicit calculations are possible for a given point set. This way it is easy to create point set optimizers using ${\displaystyle L^{2}}$ discrepancy as criteria.

It is computationally hard to find the exact value of the discrepancy of large point sets. The Erdős–Turán–Koksma inequality provides an upper bound.

Let x₁,...,x_N be points in I^s and H be an arbitrary positive integer. Then

${\displaystyle D_{N}^{*}(x_{1},\ldots ,x_{N})\leq \left({\frac {3}{2}}\right)^{s}\left({\frac {2}{H+1}}+\sum _{0<\|h\|_{\infty }\leq H}{\frac {1}{r(h)}}\left|{\frac {1}{N}}\sum _{n=1}^{N}e^{2\pi i\langle h,x_{n}\rangle }\right|\right)}$

where

${\displaystyle r(h)=\prod _{i=1}^{s}\max\{1,|h_{i}|\}\quad {\mbox{for}}\quad h=(h_{1},\ldots ,h_{s})\in \mathbb {Z} ^{s}.}$

Conjecture 1. There is a constant c_s depending only on the dimension s, such that

${\displaystyle D_{N}^{*}(x_{1},\ldots ,x_{N})\geq c_{s}{\frac {(\ln N)^{s-1}}{N}}}$

for any finite point set {x₁,...,x_N}.

Conjecture 2. There is a constant c^'_s depending only on s, such that

${\displaystyle D_{N}^{*}(x_{1},\ldots ,x_{N})\geq c'_{s}{\frac {(\ln N)^{s}}{N}}}$

for infinite number of N for any infinite sequence x₁,x₂,x₃,....

These conjectures are equivalent. They have been proved for s ≤ 2 by W. M. Schmidt. In higher dimensions, the corresponding problem is still open. The best-known lower bounds are due to Michael Lacey and collaborators.

Let s = 1. Then

${\displaystyle D_{N}^{*}(x_{1},\ldots ,x_{N})\geq {\frac {1}{2N}}}$

for any finite point set {x₁, ..., x_N}.

Let s = 2. W. M. Schmidt proved that for any finite point set {x₁, ..., x_N},

${\displaystyle D_{N}^{*}(x_{1},\ldots ,x_{N})\geq C{\frac {\log N}{N}}}$

where

${\displaystyle C=\max _{a\geq 3}{\frac {1}{16}}{\frac {a-2}{a\log a}}=0.023335\dots .}$

For arbitrary dimensions s > 1, K. F. Roth proved that

${\displaystyle D_{N}^{*}(x_{1},\ldots ,x_{N})\geq {\frac {1}{2^{4s}}}{\frac {1}{((s-1)\log 2)^{\frac {s-1}{2}}}}{\frac {\log ^{\frac {s-1}{2}}N}{N}}}$

for any finite point set {x₁, ..., x_N}. Jozef Beck [1] established a double log improvement of this result in three dimensions. This was improved by D. Bilyk and M. T. Lacey to a power of a single logarithm. The best known bound for s > 2 is due D. Bilyk and M. T. Lacey and A. Vagharshakyan.[2] For s > 2 there is a t > 0 so that

${\displaystyle D_{N}^{*}(x_{1},\ldots ,x_{N})\geq t{\frac {\log ^{{\frac {s-1}{2}}+t}N}{N}}}$

for any finite point set {x₁, ..., x_N}.

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

${\displaystyle D_{N}^{*}(x_{1},\ldots ,x_{N})\leq C{\frac {(\ln N)^{s}}{N}}.}$

where C 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 sequences, and the Sobol’ sequences. One general limitation is that construction methods can usually only guarantee the order of convergence. Practically, low discrepancy can be only achieved if N is large enough, and for large given s this minimum N can be very large. This means running a Monte-Carlo analysis with e.g. s=20 variables and N=1000 points from a low-discrepancy sequence generator may offer only a very minor accuracy improvement .

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 ${\displaystyle r_{i}}$ on ${\displaystyle [0,0.5)}$ and construct quasirandom numbers ${\displaystyle s_{i}}$ which are uniform on ${\displaystyle [0,1)}$ using:

${\displaystyle s_{i}=r_{i}}$ for ${\displaystyle i}$ odd and ${\displaystyle s_{i}=0.5+r_{i}}$ for ${\displaystyle i}$ even.

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

${\displaystyle s_{i}=s_{i-1}+0.5+r_{i}{\pmod {1}}.\,}$

That is, take the previous quasirandom number, add 0.5 and the random number, and take the result modulo 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.

Coverage of the unit square. Left for additive quasirandom numbers with c = 0.5545497..., 0.308517... Right for random numbers. From top to bottom. 10, 100, 1000, 10000 points.

Halton sequence

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₁, ..., b_s be arbitrary coprime integers greater than 1. Define

${\displaystyle x(n)=(g_{b_{1}}(n),\dots ,g_{b_{s}}(n)).}$

Then there is a constant C depending only on b₁, ..., b_s, such that sequence {x(n)}_n≥1 is a s-dimensional sequence with

${\displaystyle D_{N}^{*}(x(1),\dots ,x(N))\leq C'{\frac {(\log N)^{s}}{N}}.}$

Hammersley set

2D Hammersley set of size 256

Let b₁,...,b_s−1 be coprime positive integers greater than 1. For given s and N, the s-dimensional Hammersley set of size N is defined by[7]

${\displaystyle x(n)=\left(g_{b_{1}}(n),\dots ,g_{b_{s-1}}(n),{\frac {n}{N}}\right)}$

for n = 1, ..., N. Then

${\displaystyle D_{N}^{*}(x(1),\dots ,x(N))\leq C{\frac {(\log N)^{s-1}}{N}}}$

where C is a constant depending only on b₁, ..., b_s−1. 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 N 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 s = 2, most low-discrepancy point set generators deliver at least near-optimum discrepancies.

The points plotted below are the first 100, 1000, and 10000 elements in a sequence of the Sobol' type. For comparison, 10000 elements of a sequence of pseudorandom points are also shown. The low-discrepancy sequence was generated by TOMS algorithm 659.[9] An implementation of the algorithm in Fortran is available from Netlib.

The first 100 points in a low-discrepancy sequence of the Sobol’ type.

The first 1000 points in the same sequence. These 1000 comprise the first 100, with 900 more points.

The first 10000 points in the same sequence. These 10000 comprise the first 1000, with 9000 more points.

For comparison, here are the first 10000 points in a sequence of uniformly distributed pseudorandom numbers. Regions of higher and lower density are evident.

• Discrepancy theory
• Markov chain Monte Carlo
• Quasi-Monte Carlo method
• Sparse grid
• Systematic sampling

• Dick, Josef; Pillichshammer, Friedrich (2010). Digital Nets and Sequences: Discrepancy Theory and Quasi-Monte Carlo Integration. ISBN 978-0-521-19159-3.
• Kuipers, L.; Niederreiter, H. (2005), Uniform distribution of sequences, Dover Publications, ISBN 0-486-45019-8
• Harald Niederreiter (1992). Random Number Generation and Quasi-Monte Carlo Methods. Society for Industrial and Applied Mathematics. ISBN 0-89871-295-5.
• Drmota, Michael; Tichy, Robert F. (1997). Sequences, Discrepancies and Applications. Lecture Notes in Math. Vol. 1651. Springer. ISBN 3-540-62606-9.
• Press, William H.; Flannery, Brian P.; Teukolsky, Saul A.; Vetterling, William T. (1992). Numerical Recipes in C (2nd ed.). Cambridge University Press. see Section 7.7 for a less technical discussion of low-discrepancy sequences. ISBN 0-521-43108-5.

• Collected Algorithms of the ACM (See algorithms 647, 659, and 738.)
• Quasi-Random Sequences from the GNU Scientific Library
• Quasi-random sampling subject to constraints at FinancialMathematics.Com
• C++ generator of Sobol’ sequence
• SciPy QMC API Reference: scipy.stats.qmc