<< Tutorials Tutorials Projections >>

Low Discrepancy >> Low Discrepancy > Tutorials > Gallery

Gallery

A gallery of low discrepancy sequences.

Purpose

The goal of this page is to gather experiments showing 2D pictures of low disrepancy sequences.

Van Der Corput Sequence

The Van Der Corput base 2 sequence is just the first dimension of the Halton sequence in 1 dimension. In the following script, we generate 15 points from the Halton sequence in 1 dimension.

u=lowdisc_ldgen(15,1,"halton")

The previous script produces the following output.

      
-->u=lowdisc_ldgen(15,1,"halton")
 u  =
    0.5     
    0.25    
    0.75    
    0.125   
    0.625   
    0.375   
    0.875   
    0.0625  
    0.5625  
    0.3125  
    0.8125  
    0.1875  
    0.6875  
    0.4375  
    0.9375  
   
    

We notice that the values are organized block by block, where the size of a block is a power of two. The block #1 has x=0.5 (1 point). The block #2 has x=0.25 and x=0.75 (2 points). The block #3 has x=0.125 to x=0.875 (4 points). This is consistent with the fact that the Van Der Corput sequence that we consider here has base 2.

In order to understand how this sequence fills the interval [0,1], we analyze the elements of the sequence in more detail. The following function prints each element of the sequence, by decomposing the number n=u*2^p, where u is the low discrepancy number and p is an appropriate integer. This allows to get the binary decomposition of the number u.

function printVDC(u, p)
    nrows=size(u,"r")
    for i=1:nrows
        n=u(i)*2^p;
        d=number_tobary(n,[],[],p);
        s=strcat(string(d'),"");
        mprintf("u=%f, n=%3d, binary=0.%s\n",u(i),n,s)
    end
endfunction

The following script computes the elements of the Van Der Corput sequence, block by block: we get the first point, then two more points, then 4 points, and finally 8 points.

lds = lowdisc_new("halton");
lds = lowdisc_configure(lds,"-dimension",1);
lds = lowdisc_startup (lds);
[lds,u]=lowdisc_next(lds,1);
printVDC(u,1)
[lds,u]=lowdisc_next(lds,2);
printVDC(u,2)
[lds,u]=lowdisc_next(lds,4);
printVDC(u,3)
[lds,u]=lowdisc_next(lds,8);
printVDC(u,4)
lds = lowdisc_destroy (lds);

The previous script produces the following output. For example, u=0.5625 is decomposed in base-2 as 0.1001, since 0.5625=1/2+1/2^4.

      
u=0.500000, n=  1, binary=0.1
u=0.250000, n=  1, binary=0.01
u=0.750000, n=  3, binary=0.11
u=0.125000, n=  1, binary=0.001
u=0.625000, n=  5, binary=0.101
u=0.375000, n=  3, binary=0.011
u=0.875000, n=  7, binary=0.111
u=0.062500, n=  1, binary=0.0001
u=0.562500, n=  9, binary=0.1001
u=0.312500, n=  5, binary=0.0101
u=0.812500, n= 13, binary=0.1101
u=0.187500, n=  3, binary=0.0011
u=0.687500, n= 11, binary=0.1011
u=0.437500, n=  7, binary=0.0111
u=0.937500, n= 15, binary=0.1111
    
    

Notice that the most significant digit (the 0 or 1 just after the decimal point) is alternatively 0, 1, 0, etc... Hence, the point altenates below and over the central value x=0.5.

In order to see where the points are located in the interval [0,1], we plot the points block-by-block, where the y-coordinate corresponds to the index of the block.

scf();
lds = lowdisc_new("halton");
lds = lowdisc_configure(lds,"-dimension",1);
lds = lowdisc_startup (lds);
k=1;
imax=5;
for i=1:5
    npoints=2^(i-1);
    [lds,u]=lowdisc_next(lds,npoints);
    plot(u',i,"bo")
    xstring(u',i,string(k:k+npoints-1))
    k=k+npoints;
end
lds = lowdisc_destroy (lds);
a = gca();
a.data_bounds=[
0 0
1 imax+1
];
xtitle("Van Der Corput Sequence (base 2)",..
"X","Block Index")

The previous script produces the following output.

We see that each block fills the interval with regular spacings, but that the order of the points is so that two consecutive points in the sequence always have x=0.5 in-between.

Halton Sequence

u=lowdisc_ldgen ( 1000 , 2 , "halton" );
scf(); 
plot(u(:,1),u(:,2),"bo")
xtitle("Halton: 1000 points","X1","X2")

The previous script produces the following output.

Scrambled (RR2) Halton Sequence

The following script produces a scrambled Halton sequence, with the permutation of Kocis and Whiten (RR2).

u=lowdisc_ldgen ( 1000 , 2 , "halton-scrambled" );
scf(); 
plot(u(:,1),u(:,2),"bo")
xtitle("Scrambled Halton: 1000 points","X1","X2")

The previous script produces the following output.

Leaped Halton Sequence

The following script produces a leaped Halton sequence. The leap is computed by lowdisc_haltonsuggest.

u=lowdisc_ldgen ( 1000 , 2 , "halton-leaped" );
scf(); 
plot(u(:,1),u(:,2),"bo")
xtitle("Leaped Halton: 1000 points","X1","X2")

The previous script produces the following output.

Reverse-Halton Sequence

u=lowdisc_ldgen ( 1000 , 2 , "halton-reverse" );
scf(); 
plot(u(:,1),u(:,2),"bo")
xtitle("Reverse-Halton: 1000 points","X1","X2")

The previous script produces the following output.

Sobol Sequence

u=lowdisc_ldgen ( 1000 , 2 , "sobol" );
scf(); 
plot(u(:,1),u(:,2),"bo")
xtitle("Sobol: 1000 points","X1","X2")

The previous script produces the following output.

Faure Sequence

u=lowdisc_ldgen ( 1000 , 2 , "faure" );
scf(); 
plot(u(:,1),u(:,2),"bo")
xtitle("Faure: 1000 points","X1","X2")

The previous script produces the following output.

Niederreiter Sequence

u=lowdisc_ldgen ( 1000 , 2 , "niederreiter" );
scf(); 
plot(u(:,1),u(:,2),"bo")
xtitle("Niederreiter: 1000 points","X1","X2")

The previous script produces the following output.

Scrambled Sobol Sequence

//
// No Scrambling
//
lds=lowdisc_new("sobol");
lds=lowdisc_configure(lds,"-dimension",2);
lds = lowdisc_startup (lds);
[lds,u] = lowdisc_next (lds,100);
lds=lowdisc_destroy(lds);
//
// Owen Scrambling
//
lds=lowdisc_new("sobol");
lds=lowdisc_configure(lds,"-dimension",2);
lds=lowdisc_configure(lds,"-scrambling","Owen");
lds = lowdisc_startup (lds);
[lds,uO] = lowdisc_next (lds,100);
lds=lowdisc_destroy(lds);
//
// Faure-Tezuka Scrambling
//
lds=lowdisc_new("sobol");
lds=lowdisc_configure(lds,"-dimension",2);
lds=lowdisc_configure(lds,"-scrambling","Faure-Tezuka");
lds = lowdisc_startup (lds);
[lds,uFT] = lowdisc_next (lds,100);
lds=lowdisc_destroy(lds);
//
// Owen-Faure-Tezuka Scrambling
//
lds=lowdisc_new("sobol");
lds=lowdisc_configure(lds,"-dimension",2);
lds=lowdisc_configure(lds,"-scrambling","Owen-Faure-Tezuka");
lds = lowdisc_startup (lds);
[lds,uOFT] = lowdisc_next (lds,100);
lds=lowdisc_destroy(lds);
//
// Make a plot
//
scf();
subplot(2,2,1);
plot(u(:,1),u(:,2),"b.");
xtitle("Sobol","X1","X2");
subplot(2,2,2);
plot(uO(:,1),uO(:,2),"b.");
xtitle("Owen Scrambled Sobol","X1","X2");
subplot(2,2,3);
plot(uFT(:,1),uFT(:,2),"b.");
xtitle("Faure-Tezuka Scrambled Sobol","X1","X2");
subplot(2,2,4);
plot(uOFT(:,1),uOFT(:,2),"b.");
xtitle("Owen-Faure-Tezuka Scrambled Sobol","X1","X2");

The previous script produces the following output.

Note on 2D plots

As we have seen, in 2D the Sobol, Niederreiter (base 2) and Faure sequences are essentially the same.

i=(1:2^4-1)';
u1=lowdisc_ldgen ( 2^4-1 , 2 , "sobol" );
u2=lowdisc_ldgen ( 2^4-1 , 2 , "faure" );
u3=lowdisc_ldgen ( 2^4-1 , 2 , "niederreiter" );
disp([i,u1,u2,u3])
      
           Sobol              Faure              Niederreiter        
    n      X1        X2       X1        X2       X1        X2        
    1.     0.5       0.5      0.5       0.5      0.5       0.5       
    2.     0.75      0.25     0.25      0.75     0.25      0.75      
    3.     0.25      0.75     0.75      0.25     0.75      0.25      
    4.     0.375     0.375    0.125     0.625    0.125     0.625     
    5.     0.875     0.875    0.625     0.125    0.625     0.125     
    6.     0.625     0.125    0.375     0.375    0.375     0.375     
    7.     0.125     0.625    0.875     0.875    0.875     0.875     
    8.     0.1875    0.3125   0.0625    0.9375   0.0625    0.9375    
    9.     0.6875    0.8125   0.5625    0.4375   0.5625    0.4375    
    10.    0.9375    0.0625   0.3125    0.1875   0.3125    0.1875    
    11.    0.4375    0.5625   0.8125    0.6875   0.8125    0.6875    
    12.    0.3125    0.1875   0.1875    0.3125   0.1875    0.3125    
    13.    0.8125    0.6875   0.6875    0.8125   0.6875    0.8125    
    14.    0.5625    0.4375   0.4375    0.5625   0.4375    0.5625    
    15.    0.0625    0.9375   0.9375    0.0625   0.9375    0.0625    
   
	

The points are coming by blocks which size is a power of 2. A part of the reason is because all these sequences in two dimensions are based on the first prime, which is 2.

The Niederreiter (base 2) and Faure sequences produce exactly the same points, in the same order. The Sobol points are also the same, but come in a different order. For example, the point n=15 in Sobol sequence is the same as point n=8 in Niederreiter (base 2) sequence.

In 3D

Since the 2D Niederreiter, Sobol and Faure sequences are the same, we consider 3D sequences.

u=lowdisc_ldgen ( 1000 , 3 , "niederreiter" );
scf(); 
lowdisc_proj2d(u)
subplot(2,2,3)
xstring(0.5,0.5,"Niederreiter: 1000 points")

u=lowdisc_ldgen ( 1000 , 3 , "sobol" );
scf(); 
lowdisc_proj2d(u)
subplot(2,2,3)
xstring(0.5,0.5,"Sobol: 1000 points")

u=lowdisc_ldgen ( 1000 , 3 , "faure" );
scf(); 
lowdisc_proj2d(u)
subplot(2,2,3)
xstring(0.5,0.5,"Faure: 1000 points")

The previous script produces the following output.

We see that the projection (x1,x2) is the same for Sobol and Niederreiter (base 2), as was seen for 2D sequences. But, the projections (x1,x3) and (x2,x3) of Sobol and Niederreiter are different. Moreover, all projections of Faure sequence are now different from the other sequences. Furthermore, the projection (x1,x2) of Faure is different from the Faure 2D sequence.

Indeed, the Faure sequence uses a basis which is the smallest prime number larger than the number of dimensions. In dimension s=3, the basis used by Faure is p=3. Since Sobol and Niederreiter (base 2) use the base 2, this explains why the sequences cannot be equal anymore.


Report an issue
<< Tutorials Tutorials Projections >>