Version pdf de ce document

Low Discrepancy Sequences


March 1, 2012

Premia 14


 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

Implementation of the Low Discrepancy Sequences
QMC simulation

Quasi Monte Carlo simulation consists in approximating the integral [0,1]df(u)du by 1N- i=1Nf(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 1
N- i=1Nf(ξ i) [0,1]df(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.

1 Tore-SQRT sequences

They are d-dimensional sequences, obtained by considering the multiples of suitable irrational numbers modulo 1.
Tore sequence
It is defined by :

ξn = ({n.α }) = ({n.α1},...,{n.αd })

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

      √ ---    √ ---
α  = (  p1,...,  pd)

where (p1,,pd) are the first d prime numbers.

If α1,d are algebric, then the discrepancy satisfies:

  *        (--1--)
D n(ξ) = O  n1 -ε  ,∀ε > 0

Click there to reach the implemented part: categQMC_sqrt_implemimplementation.

2 Van der Corput and Halton sequences

2.1 Van der Corput sequence

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

        R(∑n) a
φp(n) =     --i-
         i=0 pi+1

where the coefficients ai are given by the digit expansion in base p of n :

    R∑(n)   i
n =     aip

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

        [     ]
R (n) =   log p

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

         1p log(pn)     ( logn )
D*n (ξ) ≤  ---------- = O   -----
         n  log (p)          n

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 a0 = n. Thus φp(n) = n
p 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).

2.2 Halton sequence

The Halton sequence is a d-dimensional generalization of the Van Der Corput sequence. Let (p1,,pd) be the d first prime numbers, then ξn is defined by:

ξ  = (φ  (n),...,φ  (n))
 n     p1          pd

where φpi(n) is the Van der Corput sequence in base pi.
The Halton sequence satisfies :

            d                 (   d   )
D *(ξ) ≤ 1-∏  pilog(pin) = O   log-(n)
  n      n i=1  log(pi)           n

with a constant Cd = 1d-pk--1-
This constant grows to infinity super-exponentially with dimension.
Click there to reach the implemented part: categQMC_halton_implemimplementation.

2.3 Permuted (Generalized) Halton sequence

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 pi for each one-dimensional sequence.
To break correlations between the inverse radical functions of different dimensions, we realize permutations of coefficients ai.
We consider pi, 1 i d) d permutations over {0,,pi-1} such that Πpi(0) = 0.
Each term of the permuted Halton sequence is defined by:

         Π  (a )        Π   (a   )
Spi(n ) = --pi--0--+ ⋅⋅⋅ +--piR(nR)(+n1)
           pi             pi

with n = i=0R(n)a ipi.
The global sequence is given by :

ξn = (Sp1(n),...,Spd(n))

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 p1,,pd (see Kocis and Whiten [5]).

3 Faure sequence

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, ai are the coefficients of the r-adic decomposition of n

    R (n)
     ∑     i
n =     air

We consider the following transformation T :

        R∑(n) a              R∑(n)  b
T : x =     --k--↦→  T(x ) =     --k--
        k=0 rk+1            k=0 rk+1

with bk = i=kR(n)C ika i modr and Cik denote binomial coefficients.
The coefficients bk are a permutation of the ak.
Precision ... and reference ?????

The Faure sequence is defined by using successive transformations Tk:

     (                      d-1        )
ξn =   φr(n),T (φr (n)),...,T   (φr (n ))

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

The discrepancy of the sequence satisfies :

  *       dlogd(n)
Dn (ξ) ≤ C  -------

where Cd is a constant dependent on d and r : C = -1
The constant Cd 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.

4 Generalized Faure sequence

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=0R(n)a i(n)ri.
The Generalized Faure sequence is defined by :

     ( R(n)  (1)     R(n)  (d))
ξ  = ( ∑   ξn,k-,...,∑   ξn,k)
 n     k=0 rk+1     k=0 rk+1


ξn(j,)k =     c(jk),sas(n),  j ≤ d,k ≤  R(n )

c(j) = (c k,s(j)) 0kR(n),0sR(n) and c(j) = A(j)Pj-1 where A(j) is a lower triangular inversible matrix such that (ai,l) Fr and P = (Csk) for k R(n),s R(n) is built with the binomial coefficients.

The discrepancy of the sequence satisfies:

D*n(ξ) ≤ C (d,r)log-(n)

where C(d,r) -1
2log r)d.

Click there to reach the implemented part: categQMC_gen_faure_implemimplementation.

5 Nets and (t,s)-sequences

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




6 Sobol sequence

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 = 2d. 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:

     (                                                     )
ξn =  a0V (01) ⊕ ⋅⋅⋅ ⊕ aR (n)VR(1)(n);...;a0V0(d) ⊕ ⋅⋅⋅ ⊕ aR (n)VR(d)(n)

where the V i(j) are direction numbers (expressed as binary fraction) obtained from d different primitive polynomials and ai denote the coefficients of the digit expansion of n in base b = 2, given by: n = i=0R(n)a i2i.
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:

                        (          )
           (logn )d      (logn )d+1
Dn*(ξ ) ≤ Cd--------+ O   ----------
              n              n

where Cd = --2t(d)-
d!(log2)d grows superexponentially with dimension,
and for K > 0,K-dlog-d
loglogd t(d) dlogd
 log2 + 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):

P (j) = xs(j) + b1xs(j)- 1 + ⋅⋅⋅ + bs(j)-1x + 1

and we select s(j) odd integers ci(j) such that

ci <  2i+1,   0 ≤ i < s(j)

The choice for constants ci(j) is not a easy step. Sobol’ article ????? gives some explanations about this problem.
- Once we have chosen P(j) and the ci(j) for i < s(j), we use the coefficients b i through the recurrence relation :

c(j) = 2b c(j) ⊕ 22b c(j) ⊕ 2s(j)-1b     c(j)   ⊕ 2s(j)c(j)   ⊕ c(j)
 i      1 i- 1     2 i-2          s(j)-1 i-s(j)       i-s(j)   i-s(j)

to determine the ci(j) for i s(j).
- Finally we calculate V by:

  (j)   c(ij)
V i  =  -i+1-

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

- We define a binary segment of length 2s as a set of points P i whose subscripts satisfy the inequality l2s i < (l + 1)2s where l = 0, 1,.
We divide up the s-dimensional unit cube Is by the planes x k = 1
2 into 2s multidimensional small cubes, which represent binary parallelepipeds.

- Property A: If in any binary segment of length 2s of the sequence P 0,,Pi,, 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 Is by the planes x k = 14,12,34 into 22s multidimensional small cubes. If in any binary segment of the sequence P0,,Pi, of length 22s, 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 2s (resp. 22s). In practice if s increases, it becomes difficult to verify the condition because we need to simulate at least 2s (22s) points.

Click there to reach the implemented part: categQMC_sobol_implemimplementation.

7 Niederreiter sequence

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:

     (                      )
       R∑(n)y(n1,)j     R∑(n)yn(s,)j
ξn = (     -j+1,...,    -j+1)
       j=0 b        j=0 b

with n = r=0R(n)a r(n)br and

      R (n)
y(i) =  ∑  c(i)a (n)  ∈ F
 n,j   r=0  j,r r        b

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 (cj,r(i)) is done at the beginning of the simulation.

The discrepancy of the sequence satisfies:

           ( (log n)s)
D *n(ξ) = O   --------

Construction of the cjr(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 P1,,Ps as the linear polynomials Pi(x) = x - bi where bq,,bs are distinct elements of Fb. Then the Niederreiter sequence is a (0,s)-sequence in base b and we have for 1 i s and j 1:
cjr(i) = 0 if 0 r < j - 1
cjr(i) = (r∕j - 1)b ir-j+1 if r j - 1.

Click there to reach the implemented part: categQMC_nied_implemimplementation.

8 General remarks on low discrepancy sequences


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