Introduction: Halton Sequences and Their Randomizations
The Halton sequence is a common low-discrepancy sequence used for Quasi-Monte Carlo simulations. It is based on the principle of using prime numbers as bases for each dimension, e.g., 2 is used for dimension 1, 3 for dimension 2, 5 for dimension 3, and so on. For each dimension, the indexes are converted from base 10 to the specific base associated with that dimension, the digits are then reversed, a decimal point is added before all the digits, and finally everything is converted back to base 10 to generate our samples.
Here is an example illustrating how Halton samples are generated: We have a 3 dimensional Halton object and the 29th index (30th sample) would be calculated the following way:
i = 29 = 11101_{2} = 1002_{3} = 104_{5}
\]
11101_{2} \;\to\; .10111_{2} = \tfrac{23}{32}
\]
1002_{3} \;\to\; .2001_{3} = \tfrac{55}{81}
\]
104_{5} \;\to\; .401_{5} = \tfrac{101}{125}
\]
Hence our 30th sample would be \( \left( \tfrac{23}{32}, \tfrac{55}{81}, \tfrac{101}{125} \right) \).
Like Lattice and Sobol sequences, Halton sequences too have certain randomizations available to generate the samples. Up till now, the two randomizations of Halton implemented in QMCPy were “QRNG”[1] and “OWEN”[2]. “QRNG” uses optimized fixed permutations of the digits and also adds a random digital shift to the digits while “OWEN” uses independent random permutations of the digits.
This blog discuses the implementation of three new randomizations for Halton: linear matrix scrambling (LMS), digital shift (DS), and linear matrix scrambling plus digital shift (LMS_DS). These randomizations are commonly used for digital nets, but have only recently been explored for Halton sequences. These new randomizations help provide unbiased estimates for QMC integration, and the variance under a linear matrix scramble is better than the variance when just using a random digital shift.
LMS, DS, and LMS_DS
- LMS: Linear Matrix Scrambling of Halton [3]. 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. After the indexes are converted to their base representations and a decimal point has been added before all the reversed digits, their dot product is computed with the scrambling matrix. After computing the dot product, the scrambled indexes/coefficients are converted to base 10 to generate our samples.
- DS: Digital Shift. Based on the bases, a different vector is generated for each dimension which is random between 0 and (base − 1). After the indexes are converted to their base representations and a decimal point has been added before all the reversed digits, they are added to the vector and then converted to base 10 to generate the samples.
- LMS_DS: A combination of Linear Matrix Scrambling and a Digital Shift. The Linear Matrix Scrambling of Halton is implemented but before converting to base 10, the Digital Shift is applied to the scrambled indexes/coefficients and then conversion to base 10 occurs to generate the samples.
“LMS” includes the origin as the first point in the sequence but “DS” helps prevent this, making “LMS_DS” the recommended randomized option.
Plot Examples of LMS, DS, and LMS_DS
import qmcpy as qp
dimension = 3
lms_halton = qp.Halton(dimension, randomize= 'LMS', generalize=False)
ds_halton = qp.Halton(dimension, randomize= 'DS', generalize=False)
lms_ds_halton = qp.Halton(dimension, randomize='LMS_DS', generalize=False)
fig1,ax1 = qp.plot_proj(lms_halton, n = 2**5, d_horizontal = range(dimension),
d_vertical = range(dimension), math_ind = False)
fig2,ax2 = qp.plot_proj(ds_halton, n = 2**5, d_horizontal = range(dimension),
d_vertical = range(dimension), math_ind = False)
fig3,ax3 = qp.plot_proj(lms_ds_halton, n = 2**5, d_horizontal = range(dimension),
d_vertical = range(dimension), math_ind = False)
fig1.suptitle("LMS")
fig2.suptitle("DS")
fig3.suptitle("LMS_DS")
- Figure 1: LMS
- Figure 2: DS
- Figure 3: LMS_DS
More plot examples can be seen in the Linear Matrix Scrambling and Digital Shift for Halton Notebook.
Speed Comparison Between All the Randomization Methods of Halton
import qmcpy as qp
from time import time
dimension = 25
samples = 100000
rand_options = ['QRNG', 'OWEN', 'LMS', 'DS', 'LMS_DS']
for i in range(len(rand_options)):
t_start = time()
rand_halton = qp.Halton(dimension, randomize=rand_options[i],
generalize=False).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.355226993560791 Time to generate samples for OWEN= 2.15480637550354 Time to generate samples for LMS= 1.348050832748413 Time to generate samples for DS= 1.0668678283691406 Time to generate samples for LMS_DS= 2.5354738235473633
Through this speed comparison, we can see that LMS_DS is slower than the other randomization techniques.
Conclusion
“LMS”, “DS”, and “LMS_DS” provide us with newer ways to manipulate Halton in order to get a more variety of samples. The digital shift prevents unrandomized Halton and Scrambled Halton from including origin as the point, which allows us to use these for different True Measure objects and Integral approximations.
References
- Hofert, M. & Lemieux, C. qrng: (Randomized) Quasi-Random Number Generators. R package version 0.0-7. 2019. https://CRAN.R-project.org/package=qrng.
- Owen, A. B. A randomized Halton algorithm in R. 2017. arXiv: 1706.02808 [stat.CO].
- Owen, A. B. & Pan, Z. Gain coefficients for scrambled Halton points. 2023. arXiv: 2308.08035 [math.NA].
Aadit Jain
Aadit Jain is a freshman at the University of California San Diego pursuing his B.S. in Electrical Engineering. He is really passionate about math and statistics and hopes to connect this passion with Electrical Engineering.