Low-discrepancy sequence explained

In mathematics, a low-discrepancy sequence is a sequence with the property that for all values of

, its subsequence

x_1,\ldots,x_N

has a low discrepancy.

Roughly speaking, the discrepancy of a sequence is low if the proportion of points in the sequence falling into an arbitrary set B is close to proportional to the measure of B, as would happen on average (but not for particular samples) in the case of an equidistributed sequence. Specific definitions of discrepancy differ regarding the choice of B (hyperspheres, hypercubes, etc.) and how the discrepancy for every B is computed (usually normalized) and combined (usually by taking the worst value).

Low-discrepancy sequences are also called quasirandom sequences, due to their common use as a replacement of uniformly distributed random numbers.The "quasi" modifier is used to denote more clearly that the values of a low-discrepancy sequence are neither random nor pseudorandom, but such sequences share some properties of random variables and in certain applications such as the quasi-Monte Carlo method their lower discrepancy is an important advantage.

Applications

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

Low-discrepancy sequences in numerical integration

Various methods of numerical integration can be phrased as approximating the integral of a function

in some interval, e.g. [0,1], as the average of the function evaluated at a set

\{x_1,...,x_N}\

in that interval:

	1
\int
	0

f(u)du ≈

	1
	N

	N
\sum
	i=1

f(x_i).

If the points are chosen as

x_i=i/N

, this is the rectangle rule.If the points are chosen to be randomly (or pseudorandomly) distributed, this is the Monte Carlo method.If the points are chosen as elements of a low-discrepancy sequence, this is the quasi-Monte Carlo method.A remarkable result, the Koksma–Hlawka inequality (stated below), shows that the error of such a method can be bounded by the product of two terms, one of which depends only on

, and the other one is the discrepancy of the set

\{x_1,...,x_N}\

It is convenient to construct the set

\{x_1,...,x_N}\

in such a way that if a set with

N+1

elements is constructed, the previous

elements need not be recomputed.The rectangle rule uses points set which have low discrepancy, but in general the elements must be recomputed if

is increased.Elements need not be recomputed in the random Monte Carlo method if

is increased,but the point sets do not have minimal discrepancy.By using low-discrepancy sequences we aim for low discrepancy and no need for recomputations, but actually low-discrepancy sequences can only be incrementally good on discrepancy if we allow no recomputation.

Definition of discrepancy

The discrepancy of a set

P=\{x_1,...,x_N}\

is defined, using Niederreiter's notation, as

D_N(P)=\sup_B\in\left|

	A(B;P)
	N

-λ_s(B)\right|

where

λ_s

is the

-dimensional Lebesgue measure,

A(B;P)

is the number of points in

that fall into

,and

is the set of

-dimensional intervals or boxes of the form

	s
\prod
	i=1

[a_i,b_i)=\{x\inR^s:a_i\lex_i<b_i\}

where

0\lea_i<b_i\le1

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

	s
\prod
	i=1

[0,u_i)

where

u_i

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

The two are related by

	*
D
	N

\leD_N\le2^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,

L²

-discrepancy or modified centered

L²

-discrepancy are also used intensively to compare the quality of uniform point sets. Both are much easier to calculate for large

and

The Koksma–Hlawka inequality

Let

\overline{I}^s

be the

-dimensional unit cube,

\overline{I}^s=[0,1] x … x [0,1]

.Let

have bounded variation

V(f)

\overline{I}^s

in the sense of Hardy and Krause.Then for any

x_1,\ldots,x_N

I^s=[0,1)^s=[0,1) x … x [0,1)

\left|

	1
	N

	N
\sum
	i=1

f(x_i)-

\int
	\barI^s

f(u)du\right| \leV(f)

	*
D
	N

(x_1,\ldots,x_N).

The Koksma–Hlawka inequality is sharp in the following sense: For any point set

\{x_1,\ldots,x_N\}

I^s

and any

\varepsilon>0

, there is a function

with bounded variation and

V(f)=1

such that

\left|

	1
	N

	N
\sum
	i=1

f(x_i)-

\int
	\barI^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

_1,\ldots,x_N)

The formula of Hlawka–Zaremba

Let

D=\{1,2,\ldots,d\}

. For

\emptyset ≠ u\subseteqD

wewrite

dx_u:=\prod_j\indx_j

and denote by

(x_u,1)

the point obtained from x by replacing thecoordinates not in u by

. Then

	1
	N

	N
\sum
	i=1

f(x_i)-

\int
	\barI^s

f(u)du= \sum_{\emptyset ≠}(-1)^|u|

\int
	[0,1]^\|u\|

\operatorname{disc}(x

u,1)	\partial^\|u\|
	\partialx_u

f(x_u,1)dx_u,

where

\operatorname{disc}(z)=

	1
	N

	N
\sum
	i=1

	d
\prod
	j=1

1
	[0,z_j)

(x_i,j)-

	d
\prod
	j=1

z_i

is the discrepancy function.

The L² version of the Koksma–Hlawka inequality

Applying the Cauchy–Schwarz inequality for integrals and sums to the Hlawka–Zaremba identity, we obtain an

L²

version of the Koksma–Hlawka inequality:

\left\|	1
	N

	N
\sum
	i=1

f(x_i)-

\int
	\barI^s

f(u)du\right|\le \|f\|_d\operatorname{disc}_d(\{t_i\}),

where

\operatorname{disc}_d(\{t_{i\})=\left(\sum}_{\emptyset ≠}

\int
	[0,1]^\|u\|

	2
\operatorname{disc}(x
	u,1)

	1/2
dx
	u\right)

and

\|f\|_d=\left(\sum_u\subseteq

\int		\left\|
	[0,1]^\|u\|

	\partial^\|u\|
	\partialx_u

	2
f(x
	u,1)\right\|

	1/2
dx
	u\right)

L²

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

L²

discrepancy as criteria.

The Erdős–Turán–Koksma inequality

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_1,\ldots,x_N

be points in

I^s

and

be an arbitrary positive integer. Then

	*
D
	N

(x_1,\ldots,x

N)\leq \left(	3
	2

s \left(	2
	H+1

\right)

+ \sum
	0<\\|h\\|_infty\leqH

	1	\left\|
	r(h)

	1
	N

	N
\sum
	n=1

	2\pii\langleh,x_n\rangle
e

\right| \right)

where

	smax\{1,\|h
r(h)=\prod
	i\|\} for

h=(h_1,\ldots,h

	s.

	s)\in\Z

The main conjectures

Conjecture 1. There is a constant

c_s

depending only on the dimension

, such that

	*
D
	N

(x_1,\ldots,x_N)\geq

s	(lnN)^s-1
	N

for any finite point set

{x_1,\ldots,x_N}

Conjecture 2. There is a constant

c'_s

depending only on :

, such that:

	*
D
	N

(x_1,\ldots,x_N)\geq

s	(lnN)^s
	N

for infinite number of

for any infinite sequence

x_1,x_2,x_3,\ldots

These conjectures are equivalent. They have been proved for

s\leq2

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.

Lower bounds

Let

s=1

. Then

	*(x
D
	1,\ldots,x

N)\geq	1
	2N

for any finite point set

\{x_1,...,x_N}\

Let

s=2

. W. M. Schmidt proved that for any finite point set

\{x_1,...,x_N}\

	*(x
D
	1,\ldots,x

_N)\geqC

	logN
	N

where

C=max_a\geq3

	1
	16

	a-2
	aloga

=0.023335....

For arbitrary dimensions

s>1

, K. F. Roth proved that

	*(x
D
	1,\ldots,x

N)\geq	1
	2^4s

	s-1
	2

((s-1)log2)

	s-1
	2

log

for any finite point set

\{x_1,...,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 dueD. Bilyk and M. T. Lacey and A. Vagharshakyan.^[2] For

\{x_1,...,x_N}\

there is a

t=0

so that

	*(x
D
	1,\ldots,x

_N)\geqt

	s-1	+t
	2

log

for any finite point set

\{x_1,...,x_N}\

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

	*
D
	N

(x_1,\ldots,x_N)\leqC

	(lnN)^s
	N

where

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

is large enough, and for large given s this minimum

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

r_i

[0,0.5)

and construct quasirandom numbers

s_i

which are uniform on

[0,1)

using:

s_i=r_i

for

odd and

s_i=0.5+r_i

for

even.

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

s_i=s_i-1+0.5+r_i\pmod1.

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.

Additive recurrence

For any irrational

\alpha

, the sequence

s_n=\{s₀+n\alpha\}

has discrepancy tending to

1/N

. Note that the sequence can be defined recursively by

s_n+1=(s_n+\alpha)\bmod1 .

A good value of

\alpha

gives lower discrepancy than a sequence of independent uniform random numbers.

The discrepancy can be bounded by the approximation exponent of

\alpha

. If the approximation exponent is

\mu

, then for any

\varepsilon>0

, the following bound holds:

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

By the Thue–Siegel–Roth theorem, the approximation exponent of any irrational algebraic number is 2, giving a bound of

N^{-1+\varepsilon}

above.

The recurrence relation above is similar to the recurrence relation used by a linear congruential generator, a poor-quality pseudorandom number generator:^[3]

r_i=(ar_i-1+c)\bmodm

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

with lowest discrepancy is the fractional part of the golden ratio:^[4]

	\sqrt{5
	-1}{2}

=\varphi-1 ≈ 0.618034.

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:

c=\sqrt{2}-1 ≈ 0.414214.

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:

c=\sqrt{2},\sqrt{3},\sqrt{5},\sqrt{7},\sqrt{11},\ldots

However, a set of values based on the generalised golden ratio has been shown to produce more evenly distributed points. ^[5]

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

(like

s>8

) other point set generators can offer much lower discrepancies.

van der Corput sequence

See main article: article and van der Corput sequence.

Let

	L-1
n=\sum
	k=0

	k
d
	k(n)b

be the

-ary representation of the positive integer

n\geq1

, i.e.

0\leqd_k(n)<b

. Set

g_b(n)=\sum

	L-1

	k=0

	-k-1
d
	k(n)b

Then there is a constant

depending only on

such that

(g_b(n))_n

satisfies

	*
D
	N(g

_b(1),...,g_b(N))\leqC

	logN
	N

where

	*
D
	N

is the star discrepancy.

Halton sequence

See main article: article and 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

x(n)=(g
	b₁

(n),...,g
	b_s

(n)).

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

	*
D		C'
	N(x(1),...,x(N))\leq

	(logN)^s
	N

Hammersley set

Let

b_1,\ldots,b_s-1

be coprime positive integers greater than 1. For given

and

, the

-dimensional Hammersley set of size

is defined by^[6]

x(n)=\left(g
	b₁

(n),...,g		(n),
	b_s-1

	n
	N

\right)

for

n=1,\ldots,N

. Then

	*
D		C
	N(x(1),...,x(N))\leq

	(logN)^s-1
	N

where

is a constant depending only on

b_1,\ldots,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

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.

Sobol sequence

See main article: article and Sobol sequence.

The Antonov–Saleev variant of the Sobol’ sequence generates numbers between zero and one directly as binary fractions of length

from a set of

special binary fractions,

V_i,i=1,2,...,w

called direction numbers. The bits of the Gray code of

G(i)

, are used to select direction numbers. To get the Sobol’ sequence value

s_i

take the exclusive or of the binary value of the Gray code of

with the appropriate direction number. The number of dimensions required affects the choice of

V_i

Poisson disk sampling

See main article: article. Poisson disk sampling is popular in video games to rapidly place objects in a way that appears random-lookingbut guarantees that every two points are separated by at least the specified minimum distance.^[7] 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.

Graphical examples

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.^[8] An implementation of the algorithm in Fortran is available from Netlib.

References

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

External links

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

Notes and References

A two-dimensional van Aardenne-Ehrenfest theorem in irregularities of distribution. Compositio Mathematica. 72 . 3. 1989. 269–339. 125940424. Beck. József. 1032337 . 0691.10041.
10.1016/j.jfa.2007.09.010. On the Small Ball Inequality in all dimensions . 2008 . Bilyk . Dmitriy . Lacey . Michael T. . Vagharshakyan . Armen . Journal of Functional Analysis . 254 . 9 . 2470–2502 . 14234006 . free . 0705.4619 .
Book: Knuth, Donald E. . . 2 . Chapter 3 – Random Numbers.
Web site: Malte . Skarupke. Fibonacci Hashing: The Optimization that the World Forgot. 16 June 2018 . 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.
Web site: Martin . Roberts. 2018 . The Unreasonable Effectiveness of Quasirandom Sequences.
Book: Monte Carlo Methods. Hammersley. J. M.. Handscomb. D. C.. 1964. 10.1007/978-94-009-5819-7. 978-94-009-5821-0 .
Herman Tulleken. Poisson Disk Sampling . Dev.Mag . 21 . March 2008. Herman . Tulleken . 21–25 .
10.1145/42288.214372. Algorithm 659 . 1988 . Bratley . Paul . Fox . Bennett L. . ACM Transactions on Mathematical Software . 14 . 88–100 . 17325779 . free .

Low-discrepancy sequence explained

Applications

Low-discrepancy sequences in numerical integration

Definition of discrepancy

The Koksma–Hlawka inequality

The formula of Hlawka–Zaremba

The L2 version of the Koksma–Hlawka inequality

The Erdős–Turán–Koksma inequality

The main conjectures

Lower bounds

Construction of low-discrepancy sequences

Random numbers

Additive recurrence

van der Corput sequence

Halton sequence

Hammersley set

Sobol sequence

Poisson disk sampling

Graphical examples

See also

References

External links

Notes and References

The L² version of the Koksma–Hlawka inequality