README
¶
qmc
Quasi-Monte Carlo sequences for Go: deterministic, low-discrepancy point sets that fill a unit hypercube more evenly than independent random sampling does. Sobol and Halton, each with optional randomization.
Try it in your browser → — watch a 39-dimensional Halton sequence collapse onto a diagonal ramp, and one toggle dissolve it, with the library itself compiled to WebAssembly.
No dependencies.
go get github.com/cwbudde/qmc
Usage
g, err := qmc.NewSobol(39, qmc.WithSkip(64), qmc.WithOwenScrambling(seed))
if err != nil {
return err
}
for i := 0; i < 600; i++ {
point := g.Next() // len(point) == 39, every coordinate in [0,1)
...
}
At(i) is the stateless form. It depends only on i and the generator's configuration,
never on how many points have been drawn, so a worker pool can share one sequence by
claiming indices from an atomic counter:
idx := counter.Add(1)
g.AtInto(int(idx)-1, pos)
Next, NextInto and Reset are stateful and not safe for concurrent use. At and
AtInto are.
Both generators satisfy qmc.Sequence, which carries exactly those six methods, so code
that just needs points can accept the interface and let the caller choose the sequence.
Bases() and Permutation() are Halton-only and deliberately not on it — Sobol works in
base 2 in every dimension and has neither, and the honest answer to a question that does
not apply is that you are holding the wrong type, not a zero value.
Which sequence
Sobol unless you have a reason. It is base 2 in every dimension and does not degrade as
dimensions are added, which is exactly where Halton struggles. It is capped at the 1024
dimensions the embedded direction numbers cover, unless you supply your own table with
WithDirectionNumbers.
Two things about Sobol's balance are worth knowing before you test it, because each one
makes a correct sequence look broken. The (t,m,s)-net property — the first 2^m points
landing one apiece in every elementary interval — holds on a 2^m-aligned block of raw
indices, so a stratification check wants WithSkip(2^m - 1); with the default skip of 0 all
40 of the first 40 dimensions come out unbalanced at m=8. And the D(6) direction numbers
optimise two-dimensional projections without making them all nets: of the 780 pairs among
the first 40 dimensions, 18 are balanced at every split at m=8 and 4 at m=10. Plot
dimensions 0 and 1 and you get one point per cell at every aspect ratio; plot 12 and 23 over
the same 256 points and 224 of the 256 cells of the 16x16 grid are empty while one holds
eight. Both are the correct table behaving correctly.
Halton has no dimension ceiling — primes are sieved on demand, so NewHalton(5000)
works — and its construction is simple enough to reproduce by hand, which matters if you
are migrating off an existing implementation. Above roughly twenty dimensions it has to be
randomized to be usable at all; see below.
Measured on a smooth 39-dimensional product integrand at n=4096 over ten randomization streams, against the same plain Monte Carlo baseline. Lower error is better; the multiplier is how many times more accurate than Monte Carlo on the same budget:
| generator | RMS relative error | vs Monte Carlo |
|---|---|---|
plain Monte Carlo (math/rand) |
4.3e-03 | 1x |
| Halton, random-digit scrambling | 2.4e-04 | 17.7x |
| Sobol, digital shift | 1.5e-04 | 29.5x |
| Sobol, Owen scrambling | 1.4e-04 | 32.0x |
| Halton, nested scrambling | 1.4e-04 | 31.9x |
Read that table with two caveats. It is one integrand, and a smooth product with decaying coefficients is the case nested scrambling suits best. And ten streams is not many: the same measurement over forty seeds moves random-digit Halton to 24.4x and nested Halton to 41.1x, which is a different ordering against Sobol than the ten-stream column shows. The ratios move with the integrand and with the stream count, so the test suite gates the generators at a factor of five against Monte Carlo rather than at any of these numbers.
Randomization
Each option below applies to one generator and is refused by name by the other, rather than ignored. They are mutually exclusive: a generator has one randomization or none.
| option | generator | what it does |
|---|---|---|
WithScrambling |
Halton | One digit permutation per dimension (Braaten & Weller 1979) |
WithNestedScrambling |
Halton | Uniform digit permutation per node, conditioned on the digits above it |
WithDigitalShift |
Sobol | One random word per dimension, XORed into every point |
WithOwenScrambling |
Sobol | Hash-based Owen scrambling: an independent flip at every node of the tree |
All four leave the low-discrepancy structure intact — each maps elementary intervals onto elementary intervals of the same size — and all four make the generator an RQMC sequence, so averaging over seeds yields an error estimate. Fix the seed and a run is reproducible.
Two honest caveats, both measured:
WithNestedScramblingintegrates better thanWithScramblingbut its worst-case adjacent-pair correlation over thirty seeds is 0.37 against 0.16, and it costs about 8x per point. It suits integration; it does not suit a parameter sweep, where the worst case is what you feel.WithOwenScramblingis nearly free onAtInto(370 ns/op against 360 for a digital shift at 39 dimensions) but roughly 3x onNextInto(197 against 65), because the Gray-code recurrence is precisely what cannot carry a non-linear scramble.
Web demo
examples/wasm-demo is a browser demo, published at
https://cwbudde.github.io/qmc/. Everything it shows is computed by this library compiled
to js/wasm — there is no JavaScript reimplementation of the sequence. It has two pages: a
Point Lab (scatter explorer plus digit inspector) and a Discrepancy Bench
(correlation heatmap plus convergence chart).
just run-wasm-demo
Introspection
These are Halton-only; Sobol has neither prime bases nor a permutation table.
Bases() returns the prime base of each dimension, in order — the d-th base is the d-th
prime, so Bases()[38] is 167. Permutation(d) returns the random-digit scrambling
permutation of {0…base-1} applied to dimension d, or nil for an unscrambled generator,
a nested-scrambled one (which has no table to hand out) or an out-of-range d. Both return
copies, so nothing a caller does to the returned slice can perturb the sequence.
fmt.Printf("dimension 38 uses base %d\n", g.Bases()[38])
Use scrambling above ~20 dimensions
The Halton sequence places its d-th coordinate by the radical inverse in base p_d, the
d-th prime. For a large base the first p_d points of that coordinate are simply
0, 1/p_d, 2/p_d, … — a ramp, not a sample — and two adjacent high-dimensional
coordinates ramp together.
Measured here at 39 dimensions and 600 points, which is what a parameter search over 39 knobs on a 600-evaluation budget actually asks for:
| configuration | median | p90 | worst |
|---|---|---|---|
| unscrambled, skip 64 | — | — | 0.81 |
WithScrambling, skip 64 |
0.093 | 0.13 | 0.16 |
WithNestedScrambling, skip 64 |
0.089 | 0.12 | 0.14 |
Those are absolute Pearson correlations over adjacent dimension pairs, taken over thirty seeds. Thirty rather than a handful because the statistic turns out to be high-variance: a change to the scrambling that was a pure re-instantiation, not a change of scheme, moved a five-seed worst case from 0.40 to 0.12. Quote the median and the tail, not one draw.
That statistic is also what changed WithNestedScrambling. Until recently it drew each
node's permutation from the affine family x → ax+b mod p rather than from all p!, which
is free of shuffles but leaves a tail this table used to show as a worst of 0.37 against
random-digit's 0.16. The cause was understood rather than guessed: at 600 points a
large-base coordinate varies only in its first digit, where an affine map is a ramp of
another slope rather than a scattering, so two dimensions drawing commensurate slopes ramp
together again. Drawing a genuine uniform permutation per node removes the tail, at about
forty times the cost per point and about a sixth of the integration advantage — the trade
is written out at the top of nested.go.
WithScrambling draws one uniform permutation of the digit alphabet per dimension and maps
every digit of that dimension's radical inverse through it (random-digit scrambling, Braaten
& Weller 1979). The permutations are independent across dimensions; within a dimension the
same permutation applies at every digit position. The sequence keeps its low-discrepancy
structure — a digit permutation maps each elementary interval onto another of the same size
— but the ramps are gone.
The cost is that the sequence becomes seed-dependent: this is randomized quasi-Monte Carlo
(RQMC), not plain QMC. Fix the seed and runs are reproducible again.
TestScramblingBreaksHighDimensionalCorrelation and
TestUnscrambledStillShowsTheDefect keep both halves of the table above honest.
Use a burn-in
The first Halton point is (1/2, 1/3, 1/5, 1/7, …), which sits near a corner of the box in
every coordinate with a large base. WithSkip(64) discards the first 64 points. It is
cheap and it is the standard remedy.
Dimensionality
Sobol covers 1024 dimensions from the embedded Joe & Kuo direction numbers.
WithDirectionNumbers(r io.Reader) takes a caller's own table in the same format for more —
upstream publishes the same construction out to 21201 dimensions at
https://web.maths.unsw.edu.au/~fkuo/sobol/, and new-joe-kuo-6.21201 can be passed whole.
The format and the invariants a table has to satisfy (contiguous d from 2, exactly s
direction numbers per row, every m_i odd and below 2^i, a primitive polynomial) are
documented on WithDirectionNumbers; anything failing them is refused at construction.
Direction numbers cannot be derived — the initial values come from those authors'
numerical search — so a table is the only honest option.
Halton has no fixed base table. Primes are sieved on demand, so NewHalton(500) works.
Scrambling allocates one permutation per dimension, sized by that dimension's prime, so
memory grows roughly as the sum of the first d primes, four bytes per entry. That is about
12 KB at 39 dimensions and 3.5 MB at 500, but the sum grows faster than d does: 5000
dimensions cost around 475 MB, so scrambling at that scale is a decision to make on purpose.
License
MIT.
The Sobol direction numbers in third_party/joe-kuo are Frances Y.
Kuo and Stephen Joe's, redistributed under their BSD-3 notice, which is kept alongside them.
Documentation
¶
Overview ¶
Package qmc provides quasi-Monte Carlo sequences: deterministic, low-discrepancy point sets that fill a unit hypercube more evenly than independent random sampling does.
Two sequences are implemented, both satisfying Sequence. Points are returned as coordinates in [0,1), so a caller maps them onto its own parameter ranges.
g, err := qmc.NewSobol(39, qmc.WithSkip(64), qmc.WithOwenScrambling(seed))
if err != nil {
return err
}
for i := 0; i < 600; i++ {
point := g.Next() // len(point) == 39, every coordinate in [0,1)
...
}
Which one to reach for. Sobol works in base 2 in every dimension and does not degrade as dimensions are added, which makes it the better default above a handful of dimensions; it is limited to the 1024 dimensions the embedded Joe-Kuo direction numbers cover, unless a caller supplies their own table. Halton has no dimension ceiling at all and is the one to keep if you need a sequence whose construction is simple enough to reproduce by hand, but above roughly twenty dimensions it must be randomized to be usable — its later coordinates degenerate into ramps that correlate with each other. See WithScrambling.
Every generator here is deterministic given its configuration, and that includes the randomizations: they are seeded, not sampled, so a run is reproducible across machines, architectures and Go versions. Randomizing is what makes these randomized quasi-Monte Carlo (RQMC) sequences — the estimator becomes unbiased and averaging over seeds gives an error estimate that plain QMC cannot offer.
The measured reason to use any of this, on a smooth 39-dimensional product integrand at 4096 points over ten streams: plain Monte Carlo reaches an RMS relative error of 4.3e-03, scrambled Halton 2.4e-04, Sobol with a digital shift 1.5e-04. The gap is structural — 1/n against 1/sqrt(n) convergence — and integration_test.go and sobol_integration_test.go hold it to at least a factor of five so that it stays true.
Index ¶
Examples ¶
Constants ¶
This section is empty.
Variables ¶
This section is empty.
Functions ¶
This section is empty.
Types ¶
type Halton ¶
type Halton struct {
// contains filtered or unexported fields
}
Halton generates points of the Halton sequence in a fixed number of dimensions.
A Halton generator is not safe for concurrent use through its stateful methods (Next, NextInto, Reset). At is stateless and may be called from any number of goroutines at once, which is the way to drive one shared sequence from a worker pool: have the workers claim indices from an atomic counter and call At.
func NewHalton ¶
NewHalton returns a generator over dims dimensions.
dims is bounded only by how many primes fit in memory; there is no fixed base table to run out of.
Example ¶
Draw points from a 5-dimensional sequence.
The first point of an unscrambled, unskipped Halton sequence is the classic one: 1/2, 1/3, 1/5, 1/7, 1/11 — the reciprocal of each dimension's prime base. Index 0 of the raw sequence is the all-zeros origin and is never returned, which is why counting starts here.
package main
import (
"fmt"
"github.com/cwbudde/qmc"
)
func main() {
g, err := qmc.NewHalton(5)
if err != nil {
panic(err)
}
for i := 0; i < 3; i++ {
p := g.Next()
fmt.Printf("point %d: %.4f %.4f %.4f %.4f %.4f\n", i, p[0], p[1], p[2], p[3], p[4])
}
}
Output: point 0: 0.5000 0.3333 0.2000 0.1429 0.0909 point 1: 0.2500 0.6667 0.4000 0.2857 0.1818 point 2: 0.7500 0.1111 0.6000 0.4286 0.2727
func (*Halton) At ¶
At returns point i of the sequence, counting from 0, without touching the cursor. It is the reproducible entry point: At(i) depends only on i and the generator's configuration, never on how many points have been drawn.
Point i corresponds to raw Halton index skip+1+i, so index 0 — the degenerate origin, all zeros before scrambling — is never returned.
Negative i is treated as 0.
Example ¶
At is the reproducible entry point: it depends only on the index and the generator's configuration, never on how many points have been drawn. That is what lets a worker pool claim indices from a shared counter and still reconstruct exactly which point produced which result — the property that makes a failed run re-runnable.
package main
import (
"fmt"
"github.com/cwbudde/qmc"
)
func main() {
g, err := qmc.NewHalton(3)
if err != nil {
panic(err)
}
// Draw some points through the cursor first; At must not care.
g.Next()
g.Next()
a := g.At(7)
b := g.At(7)
fmt.Printf("At(7) = %.4f %.4f %.4f\n", a[0], a[1], a[2])
fmt.Printf("again = %.4f %.4f %.4f\n", b[0], b[1], b[2])
fmt.Printf("identical = %v\n", a[0] == b[0] && a[1] == b[1] && a[2] == b[2])
}
Output: At(7) = 0.0625 0.8889 0.6400 again = 0.0625 0.8889 0.6400 identical = true
func (*Halton) AtInto ¶
AtInto is At without the allocation. As with NextInto, dst shorter than Dims() panics rather than being silently truncated.
Example ¶
AtInto writes into a caller-owned buffer and allocates nothing, which is what an optimizer's inner loop wants. The buffer must have room for Dims() coordinates; a shorter one panics rather than being silently truncated.
package main
import (
"fmt"
"github.com/cwbudde/qmc"
)
func main() {
g, err := qmc.NewHalton(4, qmc.WithSkip(64))
if err != nil {
panic(err)
}
point := make([]float64, g.Dims())
for i := 0; i < 2; i++ {
g.AtInto(i, point)
fmt.Printf("%d: %.4f %.4f %.4f %.4f\n", i, point[0], point[1], point[2], point[3])
}
}
Output: 0: 0.5078 0.7284 0.1360 0.3294 1: 0.2578 0.1728 0.3360 0.4723
func (*Halton) Bases ¶ added in v0.2.0
Bases returns the prime base of each dimension, in order: the d-th entry is the d-th prime, which is the base whose radical inverse produces coordinate d. That is the whole of what distinguishes one dimension from another, so it is also the number that explains a dimension's behaviour — dimension 38 uses base 167, and its first 167 points therefore march up a ramp in steps of 1/167 unless scrambling is on.
The returned slice is a fresh copy on every call. The generator reads its bases on every coordinate of every point, so handing out the internal slice would let a caller who sorts, truncates or edits it turn the sequence into a different one — and not visibly: the points would keep looking like plausible low-discrepancy points while no longer being the Halton sequence at all. A copy of a few hundred ints per call is not worth a failure mode nobody can see.
This exists so a UI (the WebAssembly demo in examples/wasm-demo does exactly this) can label what it is drawing without re-deriving the prime table the library already computed.
Example ¶
Bases reports the prime base behind each dimension. It is the number that explains a dimension's behaviour, which is why a UI drawing the sequence wants it: dimension 8 uses base 23, so unscrambled it needs 23 points before it stops looking like a ramp.
package main
import (
"fmt"
"github.com/cwbudde/qmc"
)
func main() {
g, err := qmc.NewHalton(10)
if err != nil {
panic(err)
}
fmt.Println(g.Bases())
}
Output: [2 3 5 7 11 13 17 19 23 29]
func (*Halton) NextInto ¶
NextInto writes the next point into dst. It allocates nothing, which matters in an optimizer's inner loop.
dst must have room for Dims() coordinates; a shorter one panics. Absorbing it instead would leave the tail coordinates holding zeros or stale values, which look like plausible positions and would steer a search silently.
func (*Halton) Permutation ¶ added in v0.2.0
Permutation returns the digit permutation applied to dimension dim, or nil when the generator is unscrambled.
With WithScrambling in effect, each dimension carries an independent uniform permutation of the digit alphabet {0..base-1} for its base, and every digit of the radical inverse — including the infinitely many leading zeros — is mapped through it. The returned slice is that permutation: entry i is the digit that digit i is rewritten to, so it has exactly Bases()[dim] entries and each of 0..base-1 appears once.
Nested scrambling (WithNestedScrambling) also returns nil, and not because it is unscrambled: it has no permutation table to hand out. Its permutations depend on the digits above the one being rewritten, so there is one per node of a tree with p^k nodes at depth k, derived on the fly and never stored — and, since nestedDigit evaluates only the entry it is asked for, most of them are never materialised in full even momentarily.
A dim outside [0, Dims()) returns nil rather than panicking, because the callers are display code walking a dimension list that may be out of step with the generator by a frame; nil is a thing a renderer can skip, a panic in that position takes the whole page down.
As with Bases, the returned slice is a fresh copy: the permutations are read on every scrambled coordinate, and a caller mutating one in place would silently corrupt every subsequent point of that dimension while the generator kept reporting the same configuration.
type Option ¶
type Option func(*settings)
An Option configures a generator at construction time. Options are applied in order and the resulting configuration is fixed for the generator's life: a sequence whose parameters could change mid-run would not be reproducible, which is the whole point of using a quasi-random sequence in the first place.
func WithDigitalShift ¶ added in v0.2.0
WithDigitalShift turns on digital shifting with the given seed: one uniform 32-bit word per dimension, XORed into every point's accumulator.
This is the cheapest randomization a digital net admits. It costs one XOR per coordinate against a word drawn at construction, which measures as 20% on AtInto at 39 dimensions and nothing at all on NextInto, where the shift is folded into the accumulator once at Reset and never touched again. Halton's digit scrambling, which has to look up a permutation for every digit of every coordinate, costs 27% on the same machine.
It buys two things. The first is an error estimate: a single QMC run gives one number with no way to say how far off it is, whereas several independent shifts give a spread that can be turned into a confidence interval. The second is the reason to use it even for a single run — a digital shift is a measure-preserving map of the unit cube onto itself that sends elementary intervals to elementary intervals, so the shifted point set is still the same (t,m,s)-net, and shifting removes the origin's special status without costing any of the structure.
What it does not do is repair a bad projection. A digital shift translates the whole net; if two dimensions' direction numbers give a poor two-dimensional projection, every shift of it is equally poor. That is what Owen scrambling is for, and why this is not the only randomization Sobol will offer.
func WithDirectionNumbers ¶ added in v0.2.0
WithDirectionNumbers supplies a Joe-Kuo direction-number table in place of the embedded one, read from r in the upstream text format: a header line, then rows of `d s a m_1 ... m_s` for d = 2, 3, 4, ...
The reason to reach for this is dimension count. The embedded table stops at 1024 because that is what fits the package's size budget; upstream publishes the same construction out to 21201, and a caller who needs more can pass the full file. It is also the way to use a different search criterion — upstream ships D(5) and D(7) alongside the D(6) set embedded here.
r goes through exactly the same parser and the same validator as the embedded table, so a file that is truncated, column-shifted or corrupted is refused at construction rather than turned into points. See validateDirectionRows for what that check does and does not prove.
The file format ¶
Upstream's files live at https://web.maths.unsw.edu.au/~fkuo/sobol/; the embedded table is the first 1024 dimensions of the D(6) family, new-joe-kuo-6.21201, and passing that file whole is the supported way to go past the embedded ceiling. The format is whitespace-separated text:
d s a m_i 2 1 0 1 3 2 1 1 3 4 3 1 1 3 1
One header line, then one row per dimension from d = 2 upward. Dimension 1 has no row: its polynomial is the empty one and all of its m_i are 1, which every Sobol implementation special-cases. A row is
d s a m_1 m_2 ... m_s
where d is the dimension, s the degree of that dimension's primitive polynomial, a the polynomial's s-1 interior coefficients packed into an integer (bit s-1-k holds the coefficient of x^(s-k)), and m_1..m_s the initial direction numbers — exactly s of them, no more and no fewer. The header is skipped if its first field is not an integer, so a hand-made file without one is accepted; refusing a file for the absence of a line nobody reads would be pedantry.
What a caller-generated table must satisfy ¶
A table that fails any of these is refused at construction, by name, and the reasoning behind each is in validateDirectionRows:
- d runs contiguously from 2, with no gaps and no repeats. A row's position in the file is what selects the dimension it is used for, so a single missing line moves every later dimension onto another dimension's polynomial — valid numbers, wrong dimension, no visible symptom.
- each row carries exactly s direction numbers.
- every m_i is odd. An even one clears the leading bit of V_i and destroys the linear independence the net property rests on.
- every m_i is below 2^i, so that m_i << (32-i) does not shift bits off the top of the word.
- the polynomial 1<<s | a<<1 | 1 is primitive over GF(2), not merely irreducible. This is the check a corrupted a field cannot pass by luck, and the one the direction-number recurrence actually depends on.
- s is between 1 and 32. Above 32 there is nowhere to put the initial direction numbers. The largest degree anywhere in the embedded 1024 dimensions is 13, so this bound only fires on a file that is not a Joe-Kuo table at all.
Nothing here proves the numbers are *good* — that they came from Joe and Kuo's search rather than from anywhere else — and nothing could. Primitivity and the m_i bounds are what makes a table usable; the quality of its two-dimensional projections is a search result, and a self-generated table that passes every check above will still have projections nobody optimised. TestDirectionTableBeyondTheEmbeddedCeiling drives a synthesised table of 1200 dimensions through this option and checks the resulting sequence is a net in every one of them, so the path above 1024 is exercised rather than assumed.
func WithNestedScrambling ¶ added in v0.2.0
WithNestedScrambling turns on nested digit scrambling with the given seed. Like WithScrambling it makes the generator a randomized QMC sequence: still low-discrepancy, but no longer identical across seeds.
It is not a free upgrade over WithScrambling and it is deliberately not the default. Measured at 39 dimensions, against random-digit scrambling:
- Integration is roughly twice as accurate. RMS relative error at n=4096 is 41x better than Monte Carlo over 40 seeds against random-digit's 24x, and 42x against 26x over 80. Over 10 seeds the measurement is too noisy to quote — it read 32x against 18x, and a variant differing only in the direction of a shuffle read 44x on the same seeds.
- Adjacent-pair correlation at small point counts is a shade better rather than worse. Over 30 seeds at 600 points the median is 0.089 against 0.093, the 90th percentile 0.123 against 0.126, and the worst 0.141 against 0.161.
- It costs about forty times as much per point. AtInto at 39 dimensions measured 20881 ns/op against 548 for random-digit scrambling and 467 unscrambled, medians of seven runs on one machine — the ratio is the part that travels. It is about 484 tree nodes per point, of which 366 are the leading-zero tails of the small bases, and each one now costs a draw from a uniform permutation rather than a table lookup.
So: reach for it when the budget is spent on an integral or an expectation, where the extra digit-level uniformity is what is being paid for, and when the integrand rather than the point count is what the wall clock is going on. Keep WithScrambling when the points are cheap to consume — a parameter sweep, a set of trial configurations — where forty times the cost per point buys an improvement in a statistic that was already acceptable.
Before this version this option drew its per-node permutations from the affine family x -> a*x+b mod p rather than from all p!. It integrated somewhat better and had a much heavier correlation tail; the points it produces have changed. See the top of nested.go for both measurements.
func WithOwenScrambling ¶ added in v0.2.0
WithOwenScrambling turns on hash-based Owen scrambling with the given seed.
It applies to Sobol generators only. Halton is not base 2 in any dimension but its first, so this construction has nothing to permute there; the nearest equivalent for Halton is WithNestedScrambling, and NewHalton says so by name rather than ignoring the option.
Prefer this to WithDigitalShift unless the cost matters. Both make the generator a randomized QMC sequence and both leave the (t,m,s)-net structure intact, but a digital shift translates the whole point set rigidly, so a poorly distributed projection stays poorly distributed under every shift it could be given. Owen scrambling actually redistributes, which is why it is the construction the theory's better convergence rates are stated for.
It subsumes the digital shift: the flip at the root of the tree is a random bit flip of the whole coordinate, which is what a one-bit digital shift is. So there is no reason to want both, and the two options are mutually exclusive rather than combinable — see the randomization type in options.go.
What it buys and what it costs, both measured at 39 dimensions. On the integrand in sobol_integration_test.go it is 1.08x more accurate than a digital shift over ten streams — a real but small margin, and small is the honest word for it on an integrand this smooth; the gap widens on functions whose projections are where the difficulty lives, which is the case Owen scrambling is for.
The cost is lopsided, and which entry point you use decides it. On AtInto it is nearly free: 369.6 ns/op against 359.9 for a digital shift, because that path already XORs one direction number per set bit of the index and a few more ALU operations disappear into it. On NextInto it is 196.6 ns/op against 65.2, a factor of three — because the Gray-code recurrence is exactly what cannot carry a non-linear scramble, so every coordinate has to be hashed on the way out and the cheap path stops being cheap. A caller drawing points with Next in an inner loop should price that before choosing.
func WithScrambling ¶
WithScrambling turns on random-digit scrambling with the given seed.
This makes the generator a randomized quasi-Monte Carlo (RQMC) sequence: still low-discrepancy, but no longer identical across seeds. That trade is deliberate. In more than roughly twenty dimensions an unscrambled Halton sequence does not fill the box at practical sample counts — its last coordinates degenerate into linear ramps that correlate with each other — and reproducibility of a sequence that is not actually filling the box is not worth much. Fix the seed and the run is reproducible again.
Example ¶
Scrambling is what makes the sequence usable above roughly twenty dimensions. It keeps the low-discrepancy structure but breaks the lockstep ramps of the high dimensions, at the cost of making the points depend on a seed — so fix the seed and the run is reproducible again.
package main
import (
"fmt"
"github.com/cwbudde/qmc"
)
func main() {
for _, seed := range []uint64{1, 2} {
g, err := qmc.NewHalton(39, qmc.WithSkip(64), qmc.WithScrambling(seed))
if err != nil {
panic(err)
}
p := g.At(0)
fmt.Printf("seed %d: dim 0 = %.4f, dim 38 = %.4f\n", seed, p[0], p[38])
}
}
Output: seed 1: dim 0 = 0.5078, dim 38 = 0.1334 seed 2: dim 0 = 0.4922, dim 38 = 0.2366
func WithSkip ¶
WithSkip discards the first n points of the underlying sequence (a burn-in).
The first few Halton points are badly placed almost by construction: point 1 sits at (1/2, 1/3, 1/5, ...), which is a corner of the box in every coordinate that has a large base. Skipping a few dozen points is the standard remedy and costs nothing.
Negative values are treated as zero.
Example ¶
WithSkip discards a burn-in. The first few Halton points sit near a corner of the box in every large-base coordinate, so a few dozen points of burn-in is the standard remedy.
package main
import (
"fmt"
"github.com/cwbudde/qmc"
)
func main() {
plain, err := qmc.NewHalton(3)
if err != nil {
panic(err)
}
skipped, err := qmc.NewHalton(3, qmc.WithSkip(64))
if err != nil {
panic(err)
}
p := plain.At(64)
s := skipped.At(0)
fmt.Printf("plain.At(64) = %.4f %.4f %.4f\n", p[0], p[1], p[2])
fmt.Printf("skipped.At(0) = %.4f %.4f %.4f\n", s[0], s[1], s[2])
}
Output: plain.At(64) = 0.5078 0.7284 0.1360 skipped.At(0) = 0.5078 0.7284 0.1360
type Sequence ¶ added in v0.2.0
type Sequence interface {
// Dims returns the number of dimensions.
Dims() int
// Next returns the next point in a freshly allocated slice.
Next() []float64
// NextInto writes the next point into dst, allocating nothing.
NextInto(dst []float64)
// Reset rewinds the cursor so the next call to Next returns point 0.
Reset()
// At returns point i, counting from 0, without touching the cursor.
At(i int) []float64
// AtInto is At without the allocation.
AtInto(i int, dst []float64)
}
A Sequence is a quasi-random point generator over a fixed number of dimensions.
This interface is deliberately small, and it deliberately excludes everything that is specific to one construction. Halton exposes Bases and Permutation; Sobol has neither, because it works in base 2 in every dimension. Putting either on this interface would force one generator to answer a question that does not apply to it, and the honest answer to an inapplicable question is not a zero value — it is that the caller is holding the wrong type. Code that needs the prime bases wants a *Halton and should say so.
What is here is the part every construction shares: how many dimensions, the stateful cursor (Next, NextInto, Reset), and the stateless index-addressed form (At, AtInto). The split matters more than it looks. At(i) depends only on i and the generator's configuration, so it is the reproducible entry point and the one safe to call concurrently; Next carries a cursor and is not.
The contract every implementation owes:
- every coordinate lies in [0,1), strictly below 1;
- At(i) is independent of how many points have been drawn;
- a dst shorter than Dims() panics rather than being truncated, because a short write leaves the tail coordinates holding stale values that look like plausible positions;
- the points depend on the configuration alone, never on GOARCH.
type Sobol ¶ added in v0.2.0
type Sobol struct {
// contains filtered or unexported fields
}
Sobol generates points of the Sobol sequence in a fixed number of dimensions, using the Joe-Kuo direction numbers.
The sequence is generated in Gray-code order, which is the order Joe and Kuo's own generator produces and the order every reference value in this package's tests was taken from. That choice is not cosmetic and it is not reversible later: point i is the direct-form point at index gray(i) = i XOR (i>>1), so Gray-code order and index order visit the same points in different sequences, and a caller who recorded outputs under one would not recognise the other. The reason to pick it is that it is the only ordering in which the stateful path can advance with a single XOR per dimension — consecutive Gray codes differ in exactly one bit, so exactly one direction number enters or leaves the accumulator. Index order would need a variable number of XORs per step and would make Next no cheaper than At, which would leave the stateful path with no reason to exist.
The reordering costs nothing that matters. Gray coding is a bijection on [0, 2^m), so the first 2^m points are the same point set either way and the (t,m,s)-net balance property — the thing the sequence is for — is untouched.
Balance holds on aligned blocks, not on any window ¶
The (t,m,s)-net property — the first 2^m points falling one apiece into every elementary interval — is a statement about a block of 2^m *raw* indices that begins on a multiple of 2^m. It is not a statement about any 2^m consecutive points, and this is the single most likely reason for a caller to conclude the sequence is broken when it is not.
Point i is raw index skip+1+i, so the default skip of 0 hands out raw indices 1..2^m, which straddles two aligned blocks and is short by exactly the one point at the far end. Measured here at 40 dimensions and m=8, over the first 256 points: with skip 0 all 40 dimensions are unbalanced, each one leaving a single interval of the 256 empty and another holding two points. It degrades from there rather than staying near-miss — with skip 100 the same run leaves up to 101 of the 256 intervals empty. With WithSkip(2^m - 1) the raw indices are 2^m..2^(m+1)-1, one aligned block, and all 40 dimensions come out exactly balanced.
So a stratification check has to choose its window: WithSkip(2^m - 1) is what this package's own balance tests construct with, and the reason is stated here rather than left as a magic constant in the tests. WithSkip chooses which aligned block you get. There is no option that makes an unaligned window a net, because no such option could exist.
Not every two-dimensional projection is a net ¶
Joe and Kuo's D(6) criterion optimises two-dimensional projections; it does not make them all t=0, and nothing about a correct table promises that it would. Measured over all 780 pairs among the first 40 dimensions, on an aligned block: 18 pairs are balanced at every split at m=8 and 4 at m=10, and (0,1) is the only pair in both lists.
The gap between a good pair and a bad one is wide enough to be alarming. Dimensions 0 and 1 put one point in every cell of every 2^a x 2^b grid with a+b = 8. Dimensions 12 and 23, over the same 256 points, leave 224 of the 256 cells of the 16x16 grid empty and pile 8 points into one cell. Both are correct output from the correct table. A caller who plots two dimensions to eyeball the sequence and happens to pick a pair like (12, 23) is looking at a real property of Sobol sequences — t grows with s, and a projection inherits no guarantee from the full-dimensional net — not at a defect. Picking a different pair is the cheap answer. A digital shift is not: it translates the whole net, so every shift of a poor projection is equally poor, which is the point WithDigitalShift's own doc comment makes.
A Sobol generator is not safe for concurrent use through its stateful methods (Next, NextInto, Reset). At and AtInto are stateless and may be called from any number of goroutines at once, which is the way to drive one shared sequence from a worker pool: have the workers claim indices from an atomic counter and call AtInto.
func NewSobol ¶ added in v0.2.0
NewSobol returns a generator over dims dimensions.
dims is limited by the direction-number table: 1024 with the embedded one, or whatever the table passed to WithDirectionNumbers covers. Exceeding it is an error rather than a wrap-around onto lower dimensions, because reusing a dimension's direction numbers would make two coordinates of every point identical — a defect that a caller integrating in a few hundred dimensions would have no way to see in the output.
func (*Sobol) At ¶ added in v0.2.0
At returns point i of the sequence, counting from 0, without touching the cursor. It is the reproducible entry point: At(i) depends only on i and the generator's configuration, never on how many points have been drawn, and it is safe to call from several goroutines at once.
Point i corresponds to raw Sobol index skip+1+i, matching Halton's convention in this package. Raw index 0 is the all-zeros origin — the same degenerate point Halton's convention exists to avoid, arrived at for a different reason: Halton's index 0 has no digits to invert, Sobol's selects no direction numbers at all. Either way it is a corner of the cube that no caller wants as their first sample, and with an unshifted generator it is exactly (0, 0, ..., 0).
The mapping from i to a raw index is what decides whether a range of points is balanced, so it is worth being explicit about here rather than only in the type doc. The (t,m,s)-net property holds over 2^m raw indices starting on a multiple of 2^m; At(0)..At(2^m-1) is that block only when skip+1 is a multiple of 2^m — which is what WithSkip(2^m - 1) arranges, and what the default skip of 0 does not. Measured at 40 dimensions and m=8, taking At(0)..At(255) with skip 0 leaves every one of the 40 dimensions with an empty interval and a doubled one. See the type doc for the rest of it; the short version is that an unaligned window of a Sobol sequence is not a net and never was.
Negative i is treated as 0.
func (*Sobol) AtInto ¶ added in v0.2.0
AtInto is At without the allocation. As with NextInto, dst shorter than Dims() panics rather than being silently truncated.
func (*Sobol) Next ¶ added in v0.2.0
Next returns the next point of the sequence in a freshly allocated slice.
func (*Sobol) NextInto ¶ added in v0.2.0
NextInto writes the next point into dst and advances the cursor. It allocates nothing.
This is the path the Gray-code ordering exists for: one XOR per dimension against a single direction number, chosen by the lowest zero bit of the counter. Compare AtInto, which has to XOR one direction number per set bit of the index — up to 32 of them.
dst must have room for Dims() coordinates; a shorter one panics. Absorbing it instead would leave the tail coordinates holding zeros or stale values, which look like plausible positions and would steer a search silently.
Source Files