March 1, 2012

1 Tore-SQRT sequences

2 Van der Corput and Halton sequences

2.1 Van der Corput sequence

2.2 Halton sequence

2.3 Permuted (Generalized) Halton sequence

3 Faure sequence

4 Generalized Faure sequence

5 Nets and (t,s)-sequences

6 Sobol sequence

7 Niederreiter sequence

8 General remarks on low discrepancy sequences

2 Van der Corput and Halton sequences

2.1 Van der Corput sequence

2.2 Halton sequence

2.3 Permuted (Generalized) Halton sequence

3 Faure sequence

4 Generalized Faure sequence

5 Nets and (t,s)-sequences

6 Sobol sequence

7 Niederreiter sequence

8 General remarks on low discrepancy sequences

Implementation of the Low Discrepancy Sequences

QMC simulation

Quasi Monte Carlo simulation consists in approximating the integral
∫
_{[0,1]d}f(u)du by ∑
_{i=1}^{N}f(u_{
i}) where {ξ_{i}} are quasi-random numbers, that means
they are generated from low-discrepancy sequences. As we already have explained
it, such sequences neither are random nor pseudo-random but deterministic and
successive values are not independent. However they satisfy good properties of
equidistribution on [0, 1]^{d} and we have that ∑
_{i=1}^{N}f(ξ_{
i}) →∫
_{[0,1]d}f(u)du.

In the following sections we describe some low discrepancy sequences. We explain
their construction and discuss some of their properties, especially on their
discrepancy.

General references about the Quasi-Monte Carlo simulation are [2], [7], [8], [6] or
[4].

The implementation of the sequences are described in the implemented
part.

They are d-dimensional sequences, obtained by considering the multiples of
suitable irrational numbers modulo 1.

∙ Tore sequence

It is defined by :

where α = (α_{1},…,α_{d}) ∈ ℝ^{d} such that (1,α_{
1},…,α_{d}) are linearly independent on ℚ.

{x} = x - [x] denotes the fractional part of x.

∙ SQRT sequence

It is a particular case of the Tore sequence with

where (p_{1},…,p_{d}) are the first d prime numbers.

If α_{1},…,α_{d} are algebric, then the discrepancy satisfies:

Click there to reach the implemented part: categQMC_sqrt_implemimplementation.

This is a one-dimensional sequence defined by the radical-inverse function φ_{p} in
base p:

where the coefficients a_{i} are given by the digit expansion in base p of n
:

R(n) denotes the maximum index for which a_{R(n)} is not equals to 0.
Its value depends on n and p by the relation p^{R(n)} ≤ n < p^{R(n)+1}, that
is:

Discrepancy of the Van der Corput sequence satisfies the following majoration:

Remarks: (ref Article of Alan Jung and Silvio Galanti)

- For n < p, there is only one positive coefficient a in the decomposition in
base p, that is a_{0} = n. Thus φ_{p}(n) = and the sequence is increasing.

- There are cycles of length p in this sequence. Each subsequence of length p
(indices kp to (k + 1)p - 1) is increasing in magnitude proportionnaly to
power of 1∕p, and covers uniformly the interval [0, 1). Consequences of this
property will be studied for multidimensional sequences (especially Halton
sequence).

The Halton sequence is a d-dimensional generalization of the Van Der Corput
sequence. Let (p_{1},…,p_{d}) be the d first prime numbers, then ξ_{n} is defined
by:

where φ_{pi}(n) is the Van der Corput sequence in base p_{i}.

The Halton sequence satisfies :

with a constant C^{d} = ∏_{
1}^{d}.

This constant grows to infinity super-exponentially with dimension.

Click there to reach the implemented part: categQMC_halton_implemimplementation.

Orthogonal projections of points from the Halton sequence show non uniform
distribution for some dimensions (see Morokoff and Caflish [6], Jung ??? or
Bratley and Fox ????). This non-uniformity is due to cycles of length p_{i} for each
one-dimensional sequence.

To break correlations between the inverse radical functions of different
dimensions, we realize permutations of coefficients a_{i}.

We consider (Π_{pi}, 1 ≤ i ≤ d) d permutations over {0,…,p_{i-1}} such that Π_{pi}(0) = 0.

Each term of the permuted Halton sequence is defined by:

with n = ∑
_{i=0}^{R(n)}a_{
i}p^{i}.

The global sequence is given by :

There is no optimal choice for the permutations. We present 3 approaches to
modify the Halton sequence

∙ An algorithm was suggested by Braaten and Weller [1] for d ≤ 16 with a
possible extension to a larger d (however with significant computation).

∙ Reverse-Radix Algorithm

An other algorithm (see Kocis and Whiten [5]) consists in reversing the binary
digits of integers, expressed using a fixed number of base 2 digits and removing
any values that are too large.

This algorithm can be applied for very large values of dimensions.

∙ Halton Sequence Leaped:

This other variant for the Halton sequence consists in using only every Lth Halton
number subject to the condition that L is a prime different from all bases p_{1},…,p_{d}
(see Kocis and Whiten [5]).

This is a d-dimensional sequence.

The Faure sequence is a permutation of the Halton sequence, but it uses the same
base r for each dimension. We choose r as the smallest odd prime integer such
that r ≥ d.

Note that the k-th dimension of a d-dimensional Faure sequence is different from
the k-th dimension of a d′-dimensional Faure sequence as soon as the base r is
different.

With usual notations, a_{i} are the coefficients of the r-adic decomposition of
n

We consider the following transformation T :

with b_{k} = ∑
_{i=k}^{R(n)}C_{
i}^{k}a_{
i} modr and C_{i}^{k} denote binomial coefficients.

The coefficients b_{k} are a permutation of the a_{k}.

Precision ... and reference ?????

The Faure sequence is defined by using successive transformations T^{k}:

where φ_{r} is the Van der Corput sequence in base r.

The discrepancy of the sequence satisfies :

where C^{d} is a constant dependent on d and r : C = ()^{d}.

The constant C^{d} tends to 0 with dimension.

The Faure sequence exhibits cycles of length r but cycles are not composed of
increasing terms, except for the first dimension. For the same dimension, the
Faure sequence has generally a smaller base than the Halton one, thus cycles are
smaller too. Because we use the smallest prime number greater than the
dimension d and not the d-th prime number.

Click there to reach the implemented part: categQMC_faure_implemimplementation.

This is a d-dimensional sequence. Let r be the smallest odd prime integer, such
that r ≥ d.

The digit expansion of n in base r is given by n = ∑
_{i=0}^{R(n)}a_{
i}(n)r^{i}.

The Generalized Faure sequence is defined by :

with

c^{(j)} = (c_{
k,s}^{(j)})_{
0≤k≤R(n),0≤s≤R(n)} and c^{(j)} = A^{(j)}P^{j-1} where A^{(j)} is a lower triangular
inversible matrix such that (a_{i,l}) ∈ F_{r} and P = (C_{s}^{k}) for k ≤ R(n),s ≤ R(n) is
built with the binomial coefficients.

The discrepancy of the sequence satisfies:

where C(d,r) ≈()^{d}.

Click there to reach the implemented part: categQMC_gen_faure_implemimplementation.

(t,s)-sequences are a group of sequences with a very regular distribution
behaviour. Their points are placed into certain equally sized volumes of the unit
cube for sequences of a fixed length. Chapter 4 of Niederreiter [7] well describes
theoretical aspects for such sequences. We just summarize in this section some
definitions and properties of those sequences.

Definitions

- An elementary interval E ∈ I
^{d}is defined as E = ∏_{ i=1}^{d}[a_{ i}b^{-di}, (a_{ i}+ 1)b^{-di}] where a_{ i},d_{i}> 0 are integers satisfying 0 ≤ a_{i}≤ b^{di}for 1 ≤ i ≤ d. - Let 0 ≤ t ≤ m be integers. A (t,m,s)-net in base b is a point set P of
b
^{m}points in I^{s}such that the number of points in E is equal to b^{t}for every elementary interval E in base b with Π(E) = b^{t-m}. - Let t ≥ 0 be an integer. A sequence x
_{0},x_{1},… of points in I^{s}is a (t,s)-sequence in base b if, for all integers k ≥ 0 and m > t, the point set constituting of the x_{n}with kb^{m}≤ n ≤ (k + 1)b^{m}is a (t,m,s)-net in base b.

Properties:

- Any (t,m,s)-net in base b is also a (u,m,s)-net in base b for integers
t ≤ u ≤ m.

The same property holds for (t,s)-sequences.

Then smaller values of t mean stronger regularity properties. - The discrepancy of a (t,m,s)-net P in base b with m > 0 satisfies:
where

- The discrepancy of the first N terms of a (t,s)-sequence P in base b
satisfies:
where

- For m ≥ 2, a (0,m,s)-net in base b can only exist if s ≤ b + 1.

A (0,s)-sequence in base b can only exists if s ≤ b.

Examples:

- The Van der Corput sequence is a (0, 1) sequence in base b. In fact, if
we consider the b
^{m}points x_{ n}with kb^{m}≤ n < (k+1)b^{m}(k ≥ 0,m ≥ 1), every b-adic interval [ab^{-m}, (a + 1)b^{-m}] contains exactly one point x_{ n}. - The s-dimensional Sobol sequence is a (τ,s)-sequence in base 2, where
τ = ∑
_{i=1}^{s}deg(P_{ i}) - s. It is called a LP_{τ}-sequence. Sobol sequence is described in the next point. - The s-dimensional Faure sequence in base r is a (0,s)-sequence where r is the smallest prime integer greater or equal than s.

The Sobol sequence is a d-dimensional sequence in base 2 and it is a (τ,d)-sequence.
It is one of the most used sequences for Quasi-Monte Carlo simulation. It was first
developped by Sobol [3] and it has been proved to have some additional
uniformity property under some initialization conditions (see [9]). Its construction
is based on primitive polynomials in the field ℤ_{2} and XOR operations.

Each dimension is a permutation of the Halton sequence with base 2
whenever N = 2^{d}. These permutations are generated from irreductible
polynomials in ℤ_{2}. But they allow for certain correlations to develop, then they
can produce regions where no points fall until N becomes very large.

The Sobol sequence is defined by:

where the V _{i}^{(j)} are direction numbers (expressed as binary fraction)
obtained from d different primitive polynomials and a_{i} denote the coefficients
of the digit expansion of n in base b = 2, given by: n = ∑
_{i=0}^{R(n)}a_{
i}2^{i}.

⊕ represents the bitwise exclusive OR operator (XOR). For explanation about
XOR operation or primitive polynomials, we refer the reader to the Numerical
Recipes in C ?????.

To implement this sequence, we use an other expression for ξ_{n} depending
only on the previous point and one direction number. This principle is
detailed in the implemented part and is due to Antonov and Saleev ?????.

The discrepancy of the sequence satisfies:

where C_{d} = grows superexponentially with dimension,

and for K > 0,K ≤ t(d) ≤ + O(d log log d). t(d) grows superlinearly
with dimension.

∙ Definition of the constants V :

- For each j ≤ d we first choose a primitive polynomial P(j) with degree
s(j):

and we select s(j) odd integers c_{i}^{(j)} such that

The choice for constants c_{i}^{(j)} is not a easy step. Sobol’ article ????? gives some
explanations about this problem.

- Once we have chosen P(j) and the c_{i}^{(j)} for i < s(j), we use the coefficients b_{
i}
through the recurrence relation :

to determine the c_{i}^{(j)} for i ≥ s(j).

- Finally we calculate V by:

∙ Uniformity property: An additional uniformity property of the sequence is
called by Sobol the property A.

- We define a binary segment of length 2^{s} as a set of points P_{
i} whose subscripts
satisfy the inequality l2^{s} ≤ i < (l + 1)2^{s} where l = 0, 1,….

We divide up the s-dimensional unit cube I^{s} by the planes x_{
k} = into 2^{s}
multidimensional small cubes, which represent binary parallelepipeds.

- Property A: If in any binary segment of length 2^{s} of the sequence P_{
0},…,P_{i},…,
all the points belong to different small cubes, then we say that the sequence
satisfies property A.

Sobol [9] proved a sufficient and necessary condition on the direction numbers so
that the property A is verified. A table of good numerical values for V is given for
a dimension s ≤ 16.

- Property A’: The property A can be extended to the property A’ defined as
follows.

We divide up the s-dimensional unit cube I^{s} by the planes x_{
k} = ,, into 2^{2s}
multidimensional small cubes. If in any binary segment of the sequence
P_{0},…,P_{i},… of length 2^{2s}, all the points belong to different small cubes, then we
say that the sequence possesses the property A’.

- Remark about the link between property A or A’ and the dimension s: Note that
the property A (resp. A’) holds for subsequences of length 2^{s} (resp. 2^{2s}). In
practice if s increases, it becomes difficult to verify the condition because we need
to simulate at least 2^{s} (2^{2s}) points.

Click there to reach the implemented part: categQMC_sobol_implemimplementation.

The Niederreiter sequence is a s-dimensional (t,s)-sequence in base b whose theoretical aspects are described in Niederreiter [7]. It is defined as:

with n = ∑
_{r=0}^{R(n)}a_{
r}(n)b^{r} and

C^{(i)} = (c_{
j,r}^{(i)}) is called the generator matrix of the i-th coordinate. An algorithm
to compute the values is given in Niederreiter [7]. Inititialization of the (c_{j,r}^{(i)}) is
done at the beginning of the simulation.

The discrepancy of the sequence satisfies:

Construction of the c_{jr}^{(i)}: (in the next version)

The method is based on the formal Laurent series.

Remark: If b is a prime power and s an arbitrary dimension such that s ≤ b, we
can choose P_{1},…,P_{s} as the linear polynomials P_{i}(x) = x - b_{i} where b_{q},…,b_{s} are
distinct elements of F_{b}. Then the Niederreiter sequence is a (0,s)-sequence in base
b and we have for 1 ≥ i ≥ s and j ≥ 1:

c_{jr}^{(i)} = 0 if 0 ≤ r < j - 1

c_{jr}^{(i)} = (r∕j - 1)b_{
i}^{r-j+1} if r ≥ j - 1.

Click there to reach the implemented part: categQMC_nied_implemimplementation.

- Quasi-random numbers combine the advantage of a random sequence that points can be added incrementally, with the advantage of a lattice that there is no clumping of points.
- For large dimension s, the theoretical bound (log N)
^{s}∕N may only be meaningful for extremly large values of N. The bound in Koksma-Hlawka inequality gives no relevant information until a very large number of points is used.

Low discrepancy sequences are very useful for low dimension. In high dimension s, a lattice can only be refined by increasing the number of points by a factor 2^{s}. - Orthogonal projections: if a d-dimensional sequence is uniformly
distributed in I
^{d}, then two-dimensional sequences formed by pairing coordinates should also be uniformly distributed. The appearance of non-uniformity in these projections is an indication of potential problems in using a quasi-random sequence for integration. This problem is developped in Morokoff and Caflish [6]. We will see that procedures like scrambling permutation can be suggested to improve the uniformity property while preserving the discrepancy.

[1] E. BRAATEN and G. WELLER. An improved low-discrepancy sequence for multidimensional quasi-monte carlo integration. Journal of Comput. Phys., (33):249–258, 1979.

[2] H.NIEDERREITER. Points sets ans sequences with small discrepancy. Monatsh.Math, 104:273–337, 1987.

[3] I.M.SOBOL. The distribution of points in a cube and the approximate evaluation of integrals. U.S.S.R. Computational Math.and Math.Phys., 7(4):86–112, 1967.

[4] INRIA. Probabilites numeriques. Chap 1: suites a discrepance faible et integration numerique.

[5] L. KOCIS and W.J. WHITEN. Computational investigations of low discrpeancy sequences. ACM Transactions on Mathematical Software, 23(2):266–294, June 1997.

[6] W.J. MOROKOFF and R.E. CAFLISH. Quasi-random sequences and their discrepancies. SIAM, Journal of Scientific Computing, 15(6):1251–1279, nov 1994.

[7] H. NIEDERREITER. Random number generation and quasi-monte carlo methods. SIAM, 1992.

[8] H. NIEDERREITER P.J.S. SHIUE. Monte carlo and quasi-monte carlo methods in scientific computing. Lecture Notes in Statistics,Ed Springer, 106, 1995.

[9] I.M. SOBOL. Uniformly distributed sequencs with an additional uniformity property. USSR Comput. Maths. Math. Phys, 16:236–242, 1976.