Low-Discrepancy Sequences and the Correlation They Leave Behind

Stratification works because it forces samples apart in a fixed set of cells. A low-discrepancy sequence generalises the idea: it produces points that are, for every prefix of the sequence, as evenly spread as any point set of that size can be. The cost of that uniformity is that the points are not independent, and the way they fail to be independent is specific and predictable. A renderer that uses such sequences without understanding the failure gets smooth noise-free images in some regions and correlated artifacts in others.

Discrepancy as a definition of even spread

The standard measure is star discrepancy: for a point set of $N$ points in the unit cube, consider every axis-aligned box anchored at the origin, and compare the fraction of points inside it with the box’s volume. The supremum of the absolute difference over all such boxes is the star discrepancy $D^_N$. A set with $D^_N = 0$ would have every anchored box exactly filled in proportion to its volume, which is impossible for a finite set; the best possible is bounded below by a constant times $\log^{d-1} N / N$ in $d$ dimensions, and sequences that achieve that bound up to a constant are called low-discrepancy.

For independent uniform random points the expected discrepancy is on the order of $N^{-1/2}$, the same as Monte Carlo error. A low-discrepancy set achieves roughly $N^{-1} \log^{d-1} N$, so for dimension 2 it is close to $N^{-1}$, and even at dimension 10 the logarithmic factor is dominated by the extra power of $N$ for any practical sample count. That is the appeal: moving from random to structured points changes the error’s exponent, not just its constant, and the improvement is large at moderate $N$.

The improvement is a worst-case guarantee, and the gap between the guarantee and observed error is where all the practical difficulties live. A function whose variation is concentrated in a small region can see error much worse than the bound suggests, and a function that is smooth sees error much better. The Koksma–Hlawka inequality makes this precise: the integration error is bounded by the function’s variation times the sequence’s discrepancy. Minimising one without the other is pointless, and the function’s variation is not something the renderer controls — it is the light transport integrand’s, which is discontinuous at every surface edge.

Building the radical inverse

The canonical construction is the van der Corput sequence in base 2. Take an integer, reverse the order of its binary digits, and write the result after a radix point. The integer 1 becomes 0.1₂ = 0.5; 2 becomes 0.01₂ = 0.25; 3 becomes 0.11₂ = 0.75. The first few elements are 1/2, 1/4, 3/4, 1/8, 5/8, 3/8, 7/8 — each new point falls into the largest remaining gap, which is the one-dimensional version of even spread.

#include <stdint.h>

/* Radical inverse in base 2: reverse the bits of n below the radix point. */
static inline double radical_inverse_2(uint32_t n)
{
    /* Bit reversal by successive shifts; O(log bits) with no lookup table. */
    n = (n << 16) | (n >> 16);
    n = ((n & 0x00FF00FFu) << 8) | ((n & 0xFF00FF00u) >> 8);
    n = ((n & 0x0F0F0F0Fu) << 4) | ((n & 0xF0F0F0F0u) >> 4);
    n = ((n & 0x33333333u) << 2) | ((n & 0xCCCCCCCCu) >> 2);
    n = ((n & 0x55555555u) << 1) | ((n & 0xAAAAAAAAu) >> 1);
    return (double)n * 2.3283064365386963e-10;   /* 2^-32 */
}

/* Halton: dimension d uses base p_d, the d-th prime. */
double halton(uint32_t n, unsigned base)
{
    double f = 1.0, r = 0.0;
    while (n > 0) {
        f /= (double)base;
        r += f * (double)(n % base);
        n /= base;
    }
    return r;
}

The bit-reversal version is exact for $2^{32}$ points and costs about ten operations per point, which is cheap enough to compute inline per sample and avoids any table or stored state. The Halton form for arbitrary bases is a loop over digits, slower but dimension-agnostic. Both are stateless functions of the sample index, which is what makes them usable in a renderer where each sample can be generated independently and in any order — a property a sequential generator does not have, and the reason these sequences are preferred for progressive rendering.

The Halton sequence’s weakness is dimension. Using the $d$-th prime as the base for dimension $d$ keeps the sequence free of correlation in the limit, but the primes grow, and in higher dimensions the first few points of a high-base sequence are concentrated at small values: the first point of base 101 is 1/101, and it takes 101 points before the sequence covers the unit interval. The early-sample structure is therefore poor in high dimensions, and the practical remedy is to skip a prefix of the sequence — a “burn-in” — or to use a different construction for high dimensions.

Why dimensions correlate

The theory of low-discrepancy sequences is entirely about points in $d$ dimensions considered jointly. That is what a renderer needs: the pixel position, the lens position for depth of field, the two direction components for a bounce, and the light sample position — all as one point in a space of perhaps 10 or 20 dimensions. The sequence provides a joint guarantee, and any subset of coordinates inherits a marginal distribution that is also well spread, but the pairs of coordinates are correlated by construction.

Here is the mechanism concretely. The van der Corput sequence’s points are the bit-reversals of consecutive integers. The bits of consecutive integers are highly structured: the least significant bit alternates, the next alternates every 2, the next every 4. Bit reversal moves the least significant bit to the most significant position of the fraction, so consecutive samples alternate between 0 and 0.5. If sample index $i$ chooses the pixel position and sample index $i$ also chooses the light position via a second dimension constructed the same way, the two are locked together: a pixel-left sample always pairs with a light-left sample. That correlation is a systematic bias, not noise, and it does not average away with more samples. It shows up as a pattern in the image: highlights that lean in one direction, shadows with an asymmetric penumbra.

The standard fix is a scramble. Applying an independent random permutation of the digits at each position, or adding a random offset modulo one in the base-$b$ sense (a digit-wise addition with carry), preserves the low-discrepancy property for each scrambled set and decorrelates the dimensions from one another. Scrambling is what makes a sequence safe to use across many dimensions: each dimension gets its own random permutation, so the joint distribution is random across dimensions while remaining structured within each. The randomised sequence is unbiased in the Monte Carlo sense, so sampling error estimates are valid, and it retains the low-discrepancy convergence within each dimension.

import random

def scrambled_radical_inverse(n, base, perms):
    """Digit-wise scramble of the radical inverse of n.
       `perms` holds one independent permutation table per digit position;
       a real Owen scramble indexes them by digit prefix as well."""
    f, r = 1.0, 0.0
    while n > 0:
        f /= base
        pos = len(perms) - 1
        while pos > 0 and f * base ** pos > 1.0:
            pos -= 1                    # which table applies to this digit
        r += f * perms[min(pos, len(perms) - 1)][n % base]
        n //= base
    return r

def build_tables(base, positions, rng):
    return [tuple(rng.sample(range(base), base)) for _ in range(positions)]

The comment is not a throwaway: a correct Owen scramble needs an independent permutation for every (position, prefix) pair, not one shared table, and the number of such permutations grows with the number of digits. Practical implementations use a hash of the prefix bits to select among a small set of tables, which trades exactness for size and is the usual compromise. The property that matters for the renderer — that the correlation between dimensions is broken — is obtained even by the cheaper variants, while the full low-discrepancy guarantee is weakened.

Structured artefacts and where they come from

Progressive refinement with an unscrambled sequence produces an image whose noise pattern changes character as samples are added. Early on, the pattern is coarse and structured by the sequence’s leading points; later it approaches uniformity. For a still image that is acceptable. For an animation with a fixed sample budget per frame, the pattern is the same in every frame, so it does not flicker as incoherent noise would; it crawls in a coherent way that is often more objectionable than noise, because the eye can track the structure. Scrambling per frame, with a different seed each frame, converts the coherent crawl into incoherent flicker, which is a real image-quality decision and not a technical detail.

The second artefact source is reusing one dimension for two purposes. A path tracer needs a random number for the choice of which BRDF lobe to sample, another for the direction within the lobe, and another for the light selection. Allocating one sequence dimension to each is correct; sharing a dimension between the bounce direction and the light position means the two are correlated in exactly the way described above, and the resulting image shows the correlation as a bias in the estimator that does not diminish with sample count. The bookkeeping is the hard part: a path tracer with variable depth needs a consistent mapping from (pixel, sample index, bounce depth, decision index) to a sequence dimension, and the mapping must not collide. Padding every bounce to a fixed number of dimensions wastes the sequence’s quality in high dimensions; using a hash of the tuple as the index into a scrambled sequence is the common working solution, and it makes the whole thing behave like stratified random sampling with the stratification of each dimension removed.

The third is sample count parity. Radiance estimates combine several samples, and if each sample contributes a different weight — as with multiple importance sampling — then the sequence’s uniformity must hold for the weighted points, not the raw ones. A low-discrepancy set that is uniform in the base measure can be badly non-uniform in the weighted measure, and the resulting estimator converges with a rate that reflects the weighted distribution, not the sequence. This is the general reason quasi-Monte Carlo does not automatically beat Monte Carlo in rendering: the integrand’s structure, the weights, and the discontinuities all have to cooperate with the sequence’s construction, and the standard practice of treating the sequence as a variance-reduction tool on top of an unbiased estimator is what keeps the results honest.

The honest comparison

The published results on quasi-Monte Carlo in rendering are consistent on one point: for low-dimensional smooth integrands — antialiasing, area light sampling with unoccluded geometry, direct illumination from a large planar emitter — structured sampling beats independent sampling decisively at moderate sample counts. For high-dimensional integrands with discontinuities — full path transport, particularly with participating media — the advantage narrows and in some configurations reverses, because the discrepancy’s logarithmic factor in dimension $d$ is not small once $d$ reaches the tens.

The practical arrangement is therefore the one that keeps the guarantee where it holds and the randomness where it protects: stratify the first few dimensions, where the structure pays, and use scrambled sequences or independent random numbers for the rest. The split point is a tuning choice, and it interacts with the convergence estimator. A quasi-Monte Carlo estimator’s error estimate cannot be computed by the sample variance, because the samples are not independent; the valid approaches are randomised replicas — run $R$ independent scrambles and take the variance across their means — or a conservative bound from the discrepancy. Treating a structured sequence’s sample variance as an error estimate is the most common misuse, and it is wrong in the direction that matters: it underestimates the error, because the sequence’s points are arranged to cancel, so the renderer reports a converged image that is still biased. The estimator’s own noise is treated in the stratified sampling article, where the variance argument has the same precondition of smoothness.

The last thing to note is sample ordering. A sequence’s low discrepancy guarantees are prefix guarantees: the first $N$ points are well distributed for every $N$. Using an arbitrary subset — the points with indices divisible by 7, say — destroys the guarantee entirely. This matters for a tiled renderer that assigns tiles to workers: each tile must either use a contiguous prefix of its own sequence or use indices computed from a global counter in a way that keeps each pixel’s set a prefix. Assigning sequence indices by tile number produces a set per tile that is a stride of the global sequence, whose discrepancy is not the sequence’s, and whose sampling pattern shows the stride as a regular offset — the very aliasing the sequence was meant to remove. The regular-grid pathology that motivates structured sampling at all is the aliasing analysed in the stratified pixel sampling article.