Add Hammersley point set discrete distribution

[Samuel-Vangu]

Motivation

QMCPy currently provides several low-discrepancy constructions, but
does not appear to provide a Hammersley point set distribution.

Hammersley point sets are a classical low-discrepancy construction
used in quasi-Monte Carlo numerical integration.

Proposed feature

Add a Hammersley point set discrete distribution following the existing
QMCPy distribution interfaces and conventions.

The construction is based on

$x_n = (\frac{n}{N}, \phi_{b_1}(n), ..., \phi_{b_{d-1}}(n)),$

where $\phi_b$ denotes the radical inverse function in base b.

Proposed implementation

• Add the new distribution to the appropriate QMCPy module.
• Follow the existing distribution API.
• Add unit tests for the mathematical construction and edge cases.
• Add documentation and examples.
• Add references to the relevant literature.

I would be interested in implementing this feature. I have previous
experience working with quasi-Monte Carlo numerical integration and
would be happy to work on it if the maintainers consider it a useful
addition.

[alegresor]

Hi @Samuel-Vangu, thanks for creating an issue! QMCPy does not currently have an explicit Hammersley class, although this would be a nice addition if you are interested in contributing!

You can get unrandomized Hammersley points in QMCPy using the Halton points and explicitly adding the first dimension, see below.

>>> import qmcpy as qp
>>> import numpy as np 
>>> d = 5 
>>> n = 8 
>>> x = np.concat([np.arange(n)[:,None]/n,qp.Halton(d-1,randomize=False,warn=False)(n)],axis=-1)
>>> x
array([[0.        , 0.        , 0.        , 0.        , 0.        ],
       [0.125     , 0.5       , 0.33333333, 0.2       , 0.14285714],
       [0.25      , 0.25      , 0.66666667, 0.4       , 0.28571429],
       [0.375     , 0.75      , 0.11111111, 0.6       , 0.42857143],
       [0.5       , 0.125     , 0.44444444, 0.8       , 0.57142857],
       [0.625     , 0.625     , 0.77777778, 0.04      , 0.71428571],
       [0.75      , 0.375     , 0.22222222, 0.24      , 0.85714286],
       [0.875     , 0.875     , 0.55555556, 0.44      , 0.02040816]])

A few potential challenges with a Hammersley implementation are

1. Extensibility: We like low-discrepancy sequences which are extensible in that increasing the number of points from the sequences does not discard previous points. This is used for example in adaptive error estimation routines which iteratively increase the number of points but do not want to discard previous queries. Hammersley is not extensible, e.g., the unrandomized sequences with n=4 and n=8 will not share the first 4 points. Some error checking for this would be required in an implementation
2. Randomization: It would be nice for a Hammersley implementation to take advantage of the existing features in our Halton implementation, especially the many randomization routines. However, I am not sure how randomization works for Hammersley or how easy this support would be to implement.

If you don't mind, please share a bit about your previous experience with QMC. Our team is always keen to connect with the QMC community and curious about both research directions and applications.

[Samuel-Vangu]

Thank you for the detailed response!

I will look into how i can handle the randomization question and the appropriate error checks.

Besides Hammersley, I am also interested in contributing other sampling constructions that could be relevant to QMCPy, including:

• Randomized Latin Hypercube Sampling;
• Niederreiter sequences;
• Sukharev grids;
• Korobov sequences.

I plan to open separate issues for these potential additions so that each one can be discussed independently.
My experience with quasi-Monte Carlo comes from a research internship at a laboratory of Université Grenoble Alpes in France, where my research topic was quasi-Monte Carlo numerical integration. During this internship, I worked with several QMC constructions and numerical integration experiments.