## QuasiMonteCarlo.jl

Lightweight and easy generation of quasi-Monte Carlo sequences with a ton of different methods on one API for easy parameter exploration in scientific machine learning (SciML)
Popularity
75 Stars
Updated Last
5 Months Ago
Started In
December 2019

# QuasiMonteCarlo.jl

This is a lightweight package for generating Quasi-Monte Carlo (QMC) samples using various different methods.

## Tutorials and Documentation

For information on using the package, see the stable documentation. Use the in-development documentation for the version of the documentation, which contains the unreleased features.

## Example

using QuasiMonteCarlo, Distributions
lb = [0.1,-0.5]
ub = [1.0,20.0]
n = 5
d = 2

s = QuasiMonteCarlo.sample(n,lb,ub,GridSample([0.1,0.5]))
s = QuasiMonteCarlo.sample(n,lb,ub,UniformSample())
s = QuasiMonteCarlo.sample(n,lb,ub,SobolSample())
s = QuasiMonteCarlo.sample(n,lb,ub,LatinHypercubeSample())
s = QuasiMonteCarlo.sample(n,lb,ub,LatticeRuleSample())
s = QuasiMonteCarlo.sample(n,lb,ub,HaltonSample([10,3], false))

The output s is a matrix, so one can use things like @uview from UnsafeArrays.jl for a stack-allocated view of the ith point:

using UnsafeArrays
@uview s[:,i]

## API

Everything has the same interface:

A = QuasiMonteCarlo.sample(n,lb,ub,sample_method)

where:

• n is the number of points to sample.
• lb is the lower bound for each variable. The length determines the dimensionality.
• ub is the upper bound.
• sample_method is the quasi-Monte Carlo sampling strategy.

Additionally, there is a helper function for generating design matrices:

k=2
As = QuasiMonteCarlo.generate_design_matrices(n,lb,ub,sample_method,k)

which returns As which is an array of k design matrices A[i] that are all sampled from the same low-discrepancy sequence.

## Available Sampling Methods

Sampling methods SamplingAlgorithm are divided into two subtypes

• DeterministicSamplingAlgorithm
• GridSample(dx) where the grid is given by lb:dx[i]:ub in the ith direction.
• SobolSample for the Sobol' sequence.
• FaureSample for the Faure sequence.
• LatticeRuleSample for a randomly-shifted rank-1 lattice rule.
• HaltonSample(base) where base[i] is the base in the ith direction.
• GoldenSample for a Golden Ratio sequence.
• KroneckerSample(alpha, s0) for a Kronecker sequence, where alpha is an length-d vector of irrational numbers (often sqrt(d)) and s0 is a length-d seed vector (often 0).
• RandomSamplingAlgorithm
• UniformSample for uniformly distributed random numbers.
• LatinHypercubeSample for a Latin Hypercube.
• Additionally, any Distribution can be used, and it will be sampled from.

## Adding a new sampling method

Adding a new sampling method is a two-step process:

1. Add a new SamplingAlgorithm type.
2. Overload the sample function with the new type.

All sampling methods are expected to return a matrix with dimension d by n, where d is the dimension of the sample space and n is the number of samples.

Example

struct NewAmazingSamplingAlgorithm{OPTIONAL} <: SamplingAlgorithm end

function sample(n,lb,ub,::NewAmazingSamplingAlgorithm)
if lb isa Number
...
return x
else
...
return reduce(hcat, x)
end
end

## Randomization of QMC sequences

Most of the previous methods are deterministic i.e. sample(n, d, Sampler()::DeterministicSamplingAlgorithm) always produces the same sequence $X = (X_1, \dots, X_n)$. The API to randomize sequence is either

• Directly use QuasiMonteCarlo.sample(n, d, DeterministicSamplingAlgorithm(R = RandomizationMethod())) or sample(n, lb, up, DeterministicSamplingAlgorithm(R = RandomizationMethod()))
• Or given any matrix $X$ ($d\times n$) of $n$ points all in dimension $d$ in $[0,1]^d$ one can do randomize(x, ::RandomizationMethod())

There are the following randomization methods:

• Scrambling methods ScramblingMethods(base, pad, rng) where base is the base used to scramble and pad the number of bits in b-ary decomposition. pad is generally chosen as $\gtrsim \log_b(n)$. The implemented ScramblingMethods are
• DigitalShift
• MatousekScramble a.k.a Linear Matrix Scramble.
• OwenScramble a.k.a Nested Uniform Scramble is the most understood theoretically but is more costly to operate.
• Shift(rng) a.k.a. Cranley Patterson Rotation.

For numerous independent randomization, use generate_design_matrices(n, d, ::DeterministicSamplingAlgorithm), ::RandomizationMethod, num_mats) where num_mats is the number of independent X generated.

### Example

Randomization of a Faure sequence with various methods.

using QuasiMonteCarlo
m = 4
d = 3
b = QuasiMonteCarlo.nextprime(d)
N = b^m # Number of points
pad = m # Length of the b-ary decomposition number = sum(y[k]/b^k for k in 1:pad)

# Unrandomized low discrepency sequence
x_faure = QuasiMonteCarlo.sample(N, d, FaureSample())

# Randomized version
x_shift = randomize(x_faure, Shift())
x_uniform = rand(d, N) # plain i.i.d. uniform
using Plots
# Setting I like for plotting
default(fontfamily="Computer Modern", linewidth=1, label=nothing, grid=true, framestyle=:default)

Plot the resulting sequences along dimensions 1 and 3.

d1 = 1
d2 = 3 # you can try every combination of dimension (d1, d2)
sequences = [x_uniform, x_faure, x_shift, x_digital_shift, x_lms, x_nus]
names = ["Uniform", "Faure (deterministic)", "Shift", "Digital Shift", "Matousek Scramble", "Owen Scramble"]
p = [plot(thickness_scaling=1.5, aspect_ratio=:equal) for i in sequences]
for (i, x) in enumerate(sequences)
scatter!(p[i], x[d1, :], x[d2, :], ms=0.9, c=1, grid=false)
title!(names[i])
xlims!(p[i], (0, 1))
ylims!(p[i], (0, 1))
yticks!(p[i], [0, 1])
xticks!(p[i], [0, 1])
hline!(p[i], range(0, 1, step=1 / 4), c=:gray, alpha=0.2)
vline!(p[i], range(0, 1, step=1 / 4), c=:gray, alpha=0.2)
hline!(p[i], range(0, 1, step=1 / 2), c=:gray, alpha=0.8)
vline!(p[i], range(0, 1, step=1 / 2), c=:gray, alpha=0.8)
end
plot(p..., size=(1500, 900)) 