Choosing the right type (Float32, Float64, etc.)

Output type

The Julia default for Floats is Float64. For QMC computations that already require a large amount of memory, this is usually not needed. In fact, like neural network packages, most QMC packages will not let you control the output type. The underlying C/C++ implementation is made with Float32 to save memory.

In this package, thanks to Julia's multiple dispatch and generic coding strategy, with a very small coding effort, one can get flexible output type, useful in different situations. For example, for research/academic purposes, it might be useful to have the exact Rational representation of a Faure or Sobol sequence, or Float16, Float64.

Example different output types

using QuasiMonteCarlo, Random
using BenchmarkTools, StatsBase
Random.seed!(1234)
m = 4
d = 2
b = QuasiMonteCarlo.nextprime(d)
N = b^m # Number of points
pad = m # Can also choose something as `2m` to get "better" randomization

QuasiMonteCarlo.sample(N, d, FaureSample(R = NoRand())) # Float64

2×16 Matrix{Float64}:
 0.015625  0.515625  0.265625  0.765625  …  0.703125  0.453125  0.953125
 0.015625  0.515625  0.765625  0.265625     0.828125  0.578125  0.078125

QuasiMonteCarlo.sample(N, d, FaureSample(R = NoRand()), Float32) # Float32

2×16 Matrix{Float32}:
 0.015625  0.515625  0.265625  0.765625  …  0.703125  0.453125  0.953125
 0.015625  0.515625  0.765625  0.265625     0.828125  0.578125  0.078125

QuasiMonteCarlo.sample(N, d, FaureSample(R = NoRand()), Rational) # Exact Rational

2×16 Matrix{Rational}:
 1//64  33//64  17//64  49//64   9//64  …  13//64  45//64  29//64  61//64
 1//64  33//64  49//64  17//64  41//64     21//64  53//64  37//64   5//64

Intermediate bits array type

For scrambling methods, large bit arrays with values in 0:base are created. Again, the default for storing is Int (i.e. generally Int64 as 64-bit Julia is the common installation) which can dramatically affect memory usage when doing multiple large computations. To avoid this, one can define the ScrambleMethod with Int32 (or else) like OwenScramble(base = Int32(b), pad = Int32(m)).^[1]

Example

To illustrate the memory gain, let's look at a bigger example with multiple randomization using iterators

m = 6
d = 3
b = QuasiMonteCarlo.nextprime(d)
N = b^m # Number of points
pad = 32 # Can also choose something as `2m` to get "better" randomization
num_mats = 50

f(x) = prod(x) * 2^length(x) # test function ∫f(x)dᵈx = 1

iterator32 = DesignMatrix(N,
    d,
    FaureSample(R = OwenScramble(base = Int32(b), pad = Int32(pad))),
    num_mats)
iterator64 = DesignMatrix(N,
    d,
    FaureSample(R = OwenScramble(base = Int(b), pad = Int(pad))),
    num_mats)
@btime [mean(f(c) for c in eachcol(X)) for X in iterator32]
# ~21.974 ms (1302 allocations: 2.18 MiB)
@btime [mean(f(c) for c in eachcol(X)) for X in iterator64]
# ~21.532 ms (1302 allocations: 3.01 MiB)

The memory usage drops from 3Mib to 2Mib! It is not exactly a factor 2 because of all the other computations happening. The number of allocations is the same, and the timing is similar. In larger cases, a time gain is also observed.