Linear Matrix Scrambling and Digital Shift for Halton¶

In Halton, each dimension has a prime number associated to it: 2 is associated to dimension 1, 3 to dimension 2, 5 to dimension 3 and so on. These prime numbers are referred as bases.

Here we explain Linear Matrix Scrambling:¶

Based on the bases, a different scrambling matrix is generated for each dimension where the lower triangle is random between 0 and (base - 1), the diagonal is random between 1 and (base - 1), and the upper triangle has all zeros. For each dimension, we convert the indices to their log base representations, compute their dot product with the scrambling matrix, and convert them to decimal (base 10) to generate our samples.

Here we explain Digital Shift:¶

Based on the bases, a different vector is generated for each dimension which is random between 0 and (base - 1). For each dimension, we convert the indices to their log base representations, add them to the vector , and convert them to decimal (base 10) to generate our samples.

Here we explain Linear Matrix Scrambling Combined with Digital Shift:¶

In this option, we implement the Linear Matrix Scrambling of Halton but before converting to base 10, we apply the Digital Shift to the scrambled coefficients (the dot product of the log base representations of the indices with the scrambling matrix), and then convert to base 10.

Here we set up the QMCPY environment:¶

In [1]:

import qmcpy as qp

import qmcpy as qp

Here we explain the parameters of the init function:¶

In [2]:

help(qp.Halton().__init__)

help(qp.Halton().__init__)

Help on method __init__ in module qmcpy.discrete_distribution.halton:

__init__(dimension=1, replications=None, seed=None, randomize='LMS_PERM', t=63, n_lim=4294967296, t_lms=None) method of qmcpy.discrete_distribution.halton.Halton instance
    Args:
        dimension (Union[int,np.ndarray]): Dimension of the generator.
    
            - If an `int` is passed in, use generating vector components at indices 0,...,`dimension`-1.
            - If an `np.ndarray` is passed in, use generating vector components at these indices.
        
        replications (int): Number of independent randomizations of a pointset.
        seed (Union[None,int,np.random.SeedSeq): Seed the random number generator for reproducibility.
        randomize (str): Options are
            
            - `'LMS_PERM'`: Linear matrix scramble with digital shift.
            - `'LMS_DS'`: Linear matrix scramble with permutation.
            - `'LMS'`: Linear matrix scramble only.
            - `'PERM'`: Permutation scramble only.
            - `'DS'`: Digital shift only.
            - `'NUS'`: Nested uniform scrambling.
            - `'QRNG'`: Deterministic permutation scramble and random digital shift from QRNG [1] (with `generalize=True`). Does *not* support replications>1.
            - `'FALSE'`: No randomization. In this case the first point will be the origin. 
        t (int): Number of bits in integer represetation of points *after* randomization. The number of bits in the generating matrices is inferred based on the largest value.
        n_lim (int): Maximum number of compatible points, determines the number of rows in the generating matrices.

Here we explain the parameters of the gen_samples function:¶

In [3]:

help(qp.Halton().gen_samples)

help(qp.Halton().gen_samples)

Help on method gen_samples in module qmcpy.discrete_distribution.abstract_discrete_distribution:

gen_samples(n=None, n_min=None, n_max=None, return_binary=False, warn=True) method of qmcpy.discrete_distribution.halton.Halton instance

Comparison between LMS, DS, and LMS_DS:¶

In [4]:

dimension = 2
lms_halton = qp.Halton(dimension, randomize= 'LMS') # Linear Matrix Scrambling
ds_halton = qp.Halton(dimension, randomize= 'DS') # Digital Shift
lms_ds_halton = qp.Halton(dimension, randomize= 'LMS_DS') # Linear Matrix Scrambling Combined with Digital Shift

dimension = 2 lms_halton = qp.Halton(dimension, randomize= 'LMS') # Linear Matrix Scrambling ds_halton = qp.Halton(dimension, randomize= 'DS') # Digital Shift lms_ds_halton = qp.Halton(dimension, randomize= 'LMS_DS') # Linear Matrix Scrambling Combined with Digital Shift

Here the difference between the three is shown by printing some samples:

In [5]:

print("Samples Generated by Linear Matrix Scrambling: ")
print(lms_halton.gen_samples(8,warn=False))
print("\nSamples Generated by Digital Shift: ")
print(ds_halton.gen_samples(8))
print("\nSamples Generated by Linear Matrix Scrambling Combined with Digital Shift: ")
print(lms_ds_halton.gen_samples(8))

print("Samples Generated by Linear Matrix Scrambling: ") print(lms_halton.gen_samples(8,warn=False)) print("\nSamples Generated by Digital Shift: ") print(ds_halton.gen_samples(8)) print("\nSamples Generated by Linear Matrix Scrambling Combined with Digital Shift: ") print(lms_ds_halton.gen_samples(8))

Samples Generated by Linear Matrix Scrambling: 
[[0.         0.        ]
 [0.66488583 0.4953392 ]
 [0.25933582 0.9904689 ]
 [0.90756417 0.21615921]
 [0.16785788 0.5632913 ]
 [0.50313275 0.72102387]
 [0.40862172 0.28368836]
 [0.76047682 0.44152192]]

Samples Generated by Digital Shift: 
[[0.06134401 0.4803077 ]
 [0.56134401 0.81364103]
 [0.31134401 0.14697436]
 [0.81134401 0.59141881]
 [0.18634401 0.92475214]
 [0.68634401 0.25808548]
 [0.43634401 0.36919659]
 [0.93634401 0.70252992]]

Samples Generated by Linear Matrix Scrambling Combined with Digital Shift: 
[[0.43966325 0.00409093]
 [0.63714101 0.72179898]
 [0.08269944 0.44092915]
 [0.77633116 0.32543015]
 [0.35612274 0.90047472]
 [0.53412944 0.60725921]
 [0.24433632 0.17094181]
 [0.92618907 0.8776755 ]]

Here the difference between the three is shown by plotting with more samples:

In [6]:

fig1,ax1 = qp.plot_proj(lms_halton, n = 2**5)
fig2,ax2 = qp.plot_proj(ds_halton,n = 2**5)
fig3,ax3 = qp.plot_proj(lms_ds_halton,n = 2**5)

fig1,ax1 = qp.plot_proj(lms_halton, n = 2**5) fig2,ax2 = qp.plot_proj(ds_halton,n = 2**5) fig3,ax3 = qp.plot_proj(lms_ds_halton,n = 2**5)

/Users/alegresor/Desktop/QMCSoftware/qmcpy/discrete_distribution/halton.py:260: UserWarning: Without randomization, the first Halton point is the origin
  warnings.warn("Without randomization, the first Halton point is the origin")

Examples of Linear Matrix Scrambling with plots:¶

Here we show a 3-dimensional scrambled Halton:

In [7]:

dimension = 3
lms_halton = qp.Halton(dimension, randomize= 'LMS')
fig,ax = qp.plot_proj(lms_halton, n = 2**5, d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False)

dimension = 3 lms_halton = qp.Halton(dimension, randomize= 'LMS') fig,ax = qp.plot_proj(lms_halton, n = 2**5, d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False)

/Users/alegresor/Desktop/QMCSoftware/qmcpy/discrete_distribution/halton.py:260: UserWarning: Without randomization, the first Halton point is the origin
  warnings.warn("Without randomization, the first Halton point is the origin")

Here we show a 4-dimensional scrambled Halton with successively increasing number of points. The initial points are in blue. The next additional points are in orange. The final additional points are in green.

In [8]:

dimension = 4
lms_halton = qp.Halton(dimension, randomize= 'LMS')
fig,ax = qp.plot_proj(lms_halton, n = [2**5,2**6,2**7], d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False,marker_size = 15)

dimension = 4 lms_halton = qp.Halton(dimension, randomize= 'LMS') fig,ax = qp.plot_proj(lms_halton, n = [2**5,2**6,2**7], d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False,marker_size = 15)

/Users/alegresor/Desktop/QMCSoftware/qmcpy/discrete_distribution/halton.py:260: UserWarning: Without randomization, the first Halton point is the origin
  warnings.warn("Without randomization, the first Halton point is the origin")

Examples of Digital Shifts with plots:¶

Here we show a 3-dimensional Halton with Digital Shift:

In [9]:

dimension = 3
ds_halton = qp.Halton(dimension, randomize= 'DS')
fig,ax = qp.plot_proj(ds_halton, n = 2**5, d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False)

dimension = 3 ds_halton = qp.Halton(dimension, randomize= 'DS') fig,ax = qp.plot_proj(ds_halton, n = 2**5, d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False)

Here we show a 4-dimensional Halton with Digital Shift with successively increasing number of points. The initial points are in blue. The next additional points are in orange. The final additional points are in green.

In [10]:

dimension = 4
ds_halton = qp.Halton(dimension, randomize= 'DS')
fig,ax = qp.plot_proj(ds_halton, n = [2**5,2**6,2**7], d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False,marker_size = 15)

dimension = 4 ds_halton = qp.Halton(dimension, randomize= 'DS') fig,ax = qp.plot_proj(ds_halton, n = [2**5,2**6,2**7], d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False,marker_size = 15)

Examples of Linear Matrix Scrambling Combined with Digital Shift with plots:¶

Here we show a 3-dimensional Halton with Linear Matrix Scrambling Combined with Digital Shift:

In [11]:

dimension = 3
lms_ds_halton = qp.Halton(dimension, randomize='LMS_DS')
fig,ax = qp.plot_proj(lms_ds_halton, n = 2**5, d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False)

dimension = 3 lms_ds_halton = qp.Halton(dimension, randomize='LMS_DS') fig,ax = qp.plot_proj(lms_ds_halton, n = 2**5, d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False)

Here we show a 4-dimensional Halton with Linear Matrix Scrambling combined with Digital Shift with successively increasing number of points. The initial points are in blue. The next additional points are in orange. The final additional points are in green.

In [12]:

dimension = 4
lms_ds_halton = qp.Halton(dimension, randomize='LMS_DS')
fig,ax = qp.plot_proj(lms_ds_halton, n = [2**5,2**6,2**7], d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False,marker_size = 15)

dimension = 4 lms_ds_halton = qp.Halton(dimension, randomize='LMS_DS') fig,ax = qp.plot_proj(lms_ds_halton, n = [2**5,2**6,2**7], d_horizontal = range(dimension), d_vertical = range(dimension),math_ind = False,marker_size = 15)

Speed Comparison Between Different Halton Randomize Options:¶

In [14]:

from time import time
dimension = 25
samples = 1000
rand_options = ['QRNG','LMS', 'DS', 'LMS_DS', 'OWEN']
for i in range(len(rand_options)):
    t_start = time()
    rand_halton = qp.Halton(dimension,randomize = rand_options[i]).gen_samples(samples,warn=False)
    t_end = time()
    print("Time to generate samples for " + rand_options[i] + " = " + str(t_end - t_start))

from time import time dimension = 25 samples = 1000 rand_options = ['QRNG','LMS', 'DS', 'LMS_DS', 'OWEN'] for i in range(len(rand_options)): t_start = time() rand_halton = qp.Halton(dimension,randomize = rand_options[i]).gen_samples(samples,warn=False) t_end = time() print("Time to generate samples for " + rand_options[i] + " = " + str(t_end - t_start))

Time to generate samples for QRNG = 0.01398015022277832
Time to generate samples for LMS = 0.02975296974182129
Time to generate samples for DS = 0.016762971878051758
Time to generate samples for LMS_DS = 0.03149604797363281
Time to generate samples for OWEN = 10.350087881088257