Low-discrepancy sequence

In mathematics, a low-discrepancy sequence is a sequence with the property that for all values of N, its subsequence x1, ..., xN 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 quasi-random or sub-random 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

Error in estimated kurtosis as a function of number of datapoints. 'Additive subrandom' 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

Subrandom numbers have an advantage over pure random numbers in that they cover the domain of interest quickly and evenly. They have an advantage over purely deterministic methods in that deterministic methods only give high accuracy when the number of datapoints is pre-set whereas in using subrandom numbers the accuracy improves continually as more datapoints are added.

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. Subrandom 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. Subrandom numbers can also be used for providing starting points for deterministic algorithms that only work locally, such as Newton–Raphson iteration.

Subrandom numbers can also be combined with search algorithms. A binary tree Quicksort-style algorithm ought to work exceptionally well because subrandom numbers flatten the tree far better than random numbers, and the flatter the tree the faster the sorting. With a search algorithm, subrandom 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

At least three methods of numerical integration can be phrased as follows. Given a set {x1, ..., xN} in the interval [0,1], approximate the integral of a function f as the average of the function evaluated at those points:

\int_0^1 f(u)\,du \approx \frac{1}{N}\,\sum_{i=1}^N f(x_i).

If the points are chosen as xi = 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 f, and the other one is the discrepancy of the set {x1, ..., xN}.

It is convenient to construct the set {x1, ..., xN} in such a way that if a set with N+1 elements is constructed, the previous N elements need not be recomputed. The rectangle rule uses points set which have low discrepancy, but in general the elements must be recomputed if N is increased. Elements need not be recomputed in the Monte Carlo method if N is increased, but the point sets do not have minimal discrepancy. By using low-discrepancy sequences, the quasi-Monte Carlo method has the desirable features of the other two methods.

Definition of discrepancy

The discrepancy of a set P = {x1, ..., xN} is defined, using Niederreiter's notation, as

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

\prod_{i=1}^s [a_i, b_i) = \{ \mathbf{x} \in \mathbf{R}^s : a_i \le x_i < b_i \} \,

where 0 \le a_i < b_i \le 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

\prod_{i=1}^s [0, u_i)

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

The two are related by

D^*_N \le D_N \le 2^s D^*_N . \,

The Koksma–Hlawka inequality

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 x1, ..., xN in Is = [0, 1) × ... × [0, 1),

\left| \frac{1}{N} \sum_{i=1}^N f(x_i) - \int_{\bar I^s} f(u)\,du \right| \le V(f)\, D_N^* (x_1,\ldots,x_N).

The Koksma-Hlawka inequality is sharp in the following sense: For any point set {x1,...,xN} in Is and any \epsilon>0, there is a function f with bounded variation and V(f)=1 such that

\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)-\epsilon.

Therefore, the quality of a numerical integration rule depends only on the discrepancy D*N(x1,...,xN).

The formula of Hlawka–Zaremba

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

dx_u:=\prod_{j\in u} dx_j

and denote by (x_u,1) the point obtained from x by replacing the coordinates not in u by 1. Then

\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|}}{\rm disc}(x_u,1)\frac{\partial^{|u|}}{\partial x_u}f(x_u,1) dx_u,

where {\rm 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.

The L^2 version of the Koksma–Hlawka inequality

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

\left|\frac{1}{N} \sum_{i=1}^N f(x_i) - \int_{\bar I^s} f(u)\,du\right|\le \|f\|_d \operatorname{disc}_d (\{t_i\}),

where

\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

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

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 x1,...,xN be points in Is and H be an arbitrary positive integer. Then

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

r(h)=\prod_{i=1}^s\max\{1,|h_i|\}\quad\mbox{for}\quad h=(h_1,\ldots,h_s)\in\Z^s.

The main conjectures

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

D_{N}^{*}(x_1,\ldots,x_N)\geq c_s\frac{(\ln N)^{s-1}}{N}

for any finite point set {x1,...,xN}.

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

D_{N}^{*}(x_1,\ldots,x_N)\geq c'_s\frac{(\ln N)^{s}}{N}

for at least one N for any infinite sequence x1,x2,x3,....

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 K. F. Roth.

Lower bounds

Let s = 1. Then

D_N^*(x_1,\ldots,x_N)\geq\frac{1}{2N}

for any finite point set {x1, ..., xN}.

Let s = 2. W. M. Schmidt proved that for any finite point set {x1, ..., xN},

D_N^*(x_1,\ldots,x_N)\geq C\frac{\log N}{N}

where

C=\max_{a\geq3}\frac{1}{16}\frac{a-2}{a\log a}=0.023335\dots

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

D_N^*(x_1,\ldots,x_N)\geq\frac{1}{2^{4s}}\frac{1}{((s-1)\log2)^\frac{s-1}{2}}\frac{\log^{\frac{s-1}{2}}N}{N}

for any finite point set {x1, ..., xN}. This bound is the best known for s > 3.

Construction of low-discrepancy sequences

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

There are constructions of sequences known such that

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.

Random numbers

Sequences of subrandom 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 on [0,0.5) and construct subrandom numbers s_i which are uniform on [0,1) using:

s_i = r_i for i odd and s_i = 0.5 + r_i for 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:

s_i = s_{i-1} + 0.5+ r_i \pmod 1. \,

That is, take the previous subrandom 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 subrandom numbers with c = 0.5545497..., 0.308517... Right for random numbers. From top to bottom. 10, 100, 1000, 10000 points.

Additive recurrence

For any irrational \alpha, the sequence

s_n = \{s_0 + n\alpha\}

has discrepancy tending to 0. (Note the sequence can be defined recursively by s_{n+1} = (s_n + \alpha)\bmod\, 1.) 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 \epsilon>0, the following bound holds:[1]

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

By the Thue-Siegel-Roth theorem, the approximation exponent of any irrational algebraic number is 2, giving a bound of N^{-1+\epsilon} above.

The value of c with lowest discrepancy is [2]

c = \frac{\sqrt{5}-1}{2} \approx 0.618034.

Another value that is nearly as good is:

c = \sqrt{2}-1 \approx 0.414214. \,

In more than one dimension, separate subrandom numbers are needed for each dimension. In higher dimensions, one set of values that can be 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 \,

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 = \bmod(a r_{i-1} + c, m)

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 list of pseudorandom number generators lists methods for generating independent pseudorandom numbers.

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 w, from a set of w special binary fractions, V_i, i = 1, 2, \dots, w called direction numbers. The bits of the Gray code of i, 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 ofi with the appropriate direction number. The number of dimensions required affects the choice of V_i.

van der Corput sequence

Let

n=\sum_{k=0}^{L-1}d_k(n)b^k

be the b-ary representation of the positive integer n ≥ 1, i.e. 0 ≤ dk(n) < b. Set

g_b(n)=\sum_{k=0}^{L-1}d_k(n)b^{-k-1}.

Then there is a constant C depending only on b such that (gb(n))n ≥ 1satisfies

D^*_N(g_b(1),\dots,g_b(N))\leq C\frac{\log N}{N},

where D*N is the star discrepancy.

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

x(n)=(g_{b_1}(n),\dots,g_{b_s}(n)).

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

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

Hammersley set

2D Hammersley set of size 256

Let b1,...,bs-1 be coprime positive integers greater than 1. For given s and N, the s-dimensional Hammersley set of size N is defined by[4]

x(n)=(g_{b_1}(n),\dots,g_{b_{s-1}}(n),\frac{n}{N})

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

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

where C is a constant depending only on b1, ..., bs−1.

Poisson disk sampling

Poisson disk sampling is popular in video games to rapidly placing objects in a way that appears random-looking but guarantees that every two points are separated by at least the specified minimum distance.[5]

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

See also


References

  1. ↑ Kuipers and Niederreiter, 2005, p. 123
  2. ↑ http://mollwollfumble.blogspot.com/, Subrandom numbers
  3. ↑ Donald E. Knuth The Art of Computer Programming Vol. 2, Ch. 3
  4. ↑ Hammersley, J. M.; Handscomb, D. C. (1964). Monte Carlo Methods. doi:10.1007/978-94-009-5819-7.
  5. ↑ Herman Tulleken. "Poisson Disk Sampling". Dev.Mag Issue 21, March 2008.
  6. ↑ P. Bratley and B.L. Fox in ACM Transactions on Mathematical Software, vol. 14, no. 1, pp 88—100

External links

This article is issued from Wikipedia - version of the Thursday, January 07, 2016. The text is available under the Creative Commons Attribution/Share Alike but additional terms may apply for the media files.