Quasi-Monte Carlo: better points,
not more points

Chapter 1: Discrepancy — measuring non-uniformity exactly

Two colleagues each hand you 100 points in the unit square and both swear theirs are "well spread." You need a single number that says which is genuinely more uniform — not a vibe, a measurement. Star discrepancy is that number.

In Chapter 0 we computed it in one dimension as a threshold scan. The full definition generalizes that to any dimension. Anchor a box at the origin (0,0) and slide its far corner to (t₁, t₂). The box has volume t₁·t₂, and some fraction of your points fall inside it. If the points were perfectly uniform, that fraction would equal the volume. The worst mismatch over all such boxes is the star discrepancy.

The word "sup" (supremum) means we hunt for the single worst test box. That is what makes star discrepancy a worst-case uniformity measure rather than an average one. A point set can look fine on average and still have one terrible box that drives its discrepancy up — and as we will see in Chapter 5, it is exactly this worst box that bounds your integration error.

Let us re-derive the two 1-D numbers from Chapter 0 as a careful warm-up, because the same scan is the heart of the canvas below. For the even midpoint set E = {1/8, 3/8, 5/8, 7/8}, the empirical step function jumps by 1/4 at each midpoint. The diagonal line y = t passes through the middle of each step, so the largest vertical gap occurs right at a jump and equals half the step height times the spacing — concretely 1/(2N) = 1/8 = 0.125.

For the clumped set C = {0.10, 0.15, 0.20, 0.70}, just below t = 0.20 we already have 3 of 4 points inside [0, t], giving a fraction of 0.75 against an interval length of 0.20 — a mismatch of 0.75 − 0.20 = 0.55. So D*(E) = 0.125 and D*(C) = 0.55. The evenly-spread set is 0.55/0.125 = 4.4× more uniform by this exact measure.

One subtlety worth internalizing: in 1-D the worst box is always achieved right at (or just after) a sample point, because the empirical fraction is a step function that only changes there. That is why our from-scratch code only needs to test the sample points, not a million thresholds. In 2-D and higher you sweep a grid of corners, which is what the auto-search button does.

python
def star_disc_2d(P, grid=200):
    N = len(P)
    ts = np.linspace(0, 1, grid+1)[1:]
    best = 0.0
    for t1 in ts:
        for t2 in ts:
            inside = np.mean((P[:,0] <= t1) & (P[:,1] <= t2))
            best = max(best, abs(inside - t1*t2))
    return best

python
from scipy.stats import qmc

# closed-form L2-star discrepancy, no grid search
qmc.discrepancy(P, method='L2-star')

Chapter 2: Van der Corput — the radical inverse

You want to drop points into [0,1] one at a time so that after any number of them, they are spread out — not just at the end. Equally-spaced points 1/N, 2/N, … cannot do this: they require you to fix N in advance and start over if you want more.

The trick is almost childish. Take the integer n, write it in a base b, then read the digits backwards behind a radix point. That single mirror operation is the radical inverse φ_b(n), and it produces a sequence that keeps bisecting every gap it has left.

Formally, write n = a₀ + a₁b + a₂b² + … in base b. The radical inverse reflects those digits across the point: φ_b(n) = a₀/b + a₁/b² + a₂/b³ + …. The least-significant digit a₀ (which alternates fastest as n increments) becomes the most-significant fractional digit, so consecutive n land on opposite ends of the interval.

Let us compute several points by hand in base 2. Take n = 6. In binary 6 = 110, so the digits are a₀ = 0, a₁ = 1, a₂ = 1. Reflect them: φ₂(6) = 0/2 + 1/4 + 1/8 = 0 + 0.25 + 0.125 = 0.375 = 3/8.

Now n = 5. In binary 5 = 101, so a₀ = 1, a₁ = 0, a₂ = 1. Reflect: φ₂(5) = 1/2 + 0/4 + 1/8 = 0.5 + 0 + 0.125 = 0.625 = 5/8. Notice that 5 and 6 land on opposite halves of [0,1] — 0.625 is in the upper half, 0.375 in the lower — even though the integers differ by only 1. That is the digit-reversal at work: it scatters neighbors.

The trick is not tied to base 2. In base 3, take n = 7. Write 7 in base 3: 7 = 2·3 + 1, so the digits are a₀ = 1, a₁ = 2 (i.e. 7 = "21" in base 3). Reflect: φ₃(7) = 1/3 + 2/9 = 3/9 + 2/9 = 5/9 = 0.5556. We will use base-3 radical inverses for the second coordinate of the Halton sequence next chapter.

Why does this fill gaps? Watch the base-2 sequence in order: φ₂(1) = 0.5, then 0.25, then 0.75, then 0.125, 0.625, 0.375, 0.875, …. The first point halves the interval. The next two land in the centers of the two halves. The next four land in the centers of the four quarters. Each new point falls in the largest current gap and bisects it — equidistribution that refines automatically, with no N fixed in advance.

python
def radical_inverse(n, base):
    f, r = 1.0, 0.0
    while n > 0:
        f /= base
        r += f * (n % base)
        n //= base
    return r

# van der Corput sequence, base 2
vdc = [radical_inverse(n, 2) for n in range(1, 9)]
# -> 0.5, 0.25, 0.75, 0.125, 0.625, 0.375, 0.875, ...

python
from scipy.stats import qmc

# Sobol's 1-D coords ARE the base-2
# van der Corput points
qmc.Sobol(d=1, scramble=False).random(8)

Chapter 3: Halton — van der Corput in many dimensions

You are placing sample points for a renderer's pixel jitter and depth-of-field — a 2-D problem — and you want them to tile the square smoothly, not in stripes. You have the van der Corput sequence from Chapter 2, but it only fills one dimension. What do you use for the second coordinate?

The naive answer fails badly. If you used base 2 for both coordinates, then point n would be (φ₂(n), φ₂(n)) — identical coordinates — and every point would sit on the diagonal line x = y. You would cover a line, not a square.

The Halton sequence fixes this by giving each axis its own prime base. Coordinate 1 uses base 2, coordinate 2 uses base 3, coordinate 3 uses base 5, and so on through the primes. Distinct (coprime) bases keep the per-dimension radical inverses decorrelated, so the joint point set fills the whole hypercube instead of marching in lock-step.

Let us compute Halton point n = 4 with bases (2, 3) by hand. Coordinate 1 is φ₂(4): 4 = 100 in binary, so digits a₀ = 0, a₁ = 0, a₂ = 1; reflect to get 0/2 + 0/4 + 1/8 = 0.125. Coordinate 2 is φ₃(4): 4 = "11" in base 3, so a₀ = 1, a₁ = 1; reflect to get 1/3 + 1/9 = 3/9 + 1/9 = 4/9 = 0.4444. So Halton point #4 is (0.125, 0.4444).

That decorrelation is not just a hope — you can measure it. For the first 16 Halton points with bases (2,3), the Pearson correlation between the two coordinate columns is about −0.19: essentially uncorrelated, nicely spread. The points fill the square.

But Halton has a famous weakness. With large prime bases, the radical inverse cycles slowly, so for the first many points two large-base coordinates move almost in step. With bases (17, 19), the correlation of the first 16 points is +1.00 — the points are essentially collinear, a thin diagonal streak before they slowly spread. This is Halton degradation, and it is why Halton is rarely used beyond about 6 dimensions without scrambling.

The fix is either to scramble the digits (which breaks the lock-step) or to switch to the Sobol sequence (next chapter), which works entirely in base 2 and stays well-distributed in much higher dimensions. The canvas lets you switch bases and literally watch the diagonal streak appear with (17, 19).

python
def halton(N, dims_bases):
    def phi(n, b):
        f, r = 1.0, 0.0
        while n > 0:
            f /= b; r += f*(n % b); n //= b
        return r
    return np.array([[phi(n, b) for b in dims_bases]
                     for n in range(1, N+1)])

P = halton(64, [2, 3])   # 64 points in [0,1]^2

python
from scipy.stats import qmc

# scrambled Halton avoids large-prime degradation
qmc.Halton(d=2, scramble=True).random(64)

Chapter 4: Sobol — digital nets done right

You need 4096 sample points in 20 dimensions for a path tracer, and Halton's high-prime axes have gone stripey. Sobol is the industry answer. Instead of mixing primes, it works entirely in base 2 with carefully chosen direction numbers, so that every dyadic box of a fixed size contains exactly its fair share of points.

That "fair share" property has a precise name: a (t, m, d)-net. Split the d-dimensional cube into elementary dyadic intervals — boxes whose side lengths are powers of 1/2 and whose volume is 2^t−m. A (t, m, d)-net guarantees that every such box contains exactly 2^t of the 2^m points. The smaller the t, the stronger the equidistribution; the best nets have t = 0.

This is the structural upgrade over Halton. Halton hopes that coprime bases keep coordinates decorrelated, and that hope fails for large primes. Sobol guarantees balance in every elementary box across all dimensions, by construction, with no primes to run out of.

Let us verify the net property by hand on scipy's first four 2-D Sobol points. The unscrambled 1-D Sobol points (the same base-2 van der Corput values, in Gray-code order) are 0, 0.5, 0.75, 0.25, 0.375, 0.875, 0.625, 0.125 — notice how 0.75 and 0.25 fill the two halves the instant they appear. The first four 2-D points are (0,0), (0.5,0.5), (0.75,0.25), (0.25,0.75).

Now bin those four points into the 2×2 dyadic grid (each cell has area 1/4). Point (0,0) lands in cell (0,0); (0.5,0.5) lands in cell (1,1); (0.75,0.25) lands in cell (1,0); (0.25,0.75) lands in cell (0,1). Each of the four cells gets exactly one point — the net property holds exactly. That is t = 0: perfect stratification of four points across four boxes.

The defining trick is that this balance holds for every elementary grid of the same total volume at once — the 1×4, the 2×2, and the 4×1 partitions all get one point per box. Random points cannot promise that: a random quadruple typically leaves a box empty and crowds another with 2 or 3. The canvas lets you flip between the three grids and confirm the green check on each.

python
# base 2, Gray-code recurrence
def sobol_1d(N, direction):  # direction = v_i ints
    x, pts, BITS = 0, [], 30
    for i in range(1, N+1):
        c = (i ^ (i-1)).bit_length() - 1  # low zero bit
        x ^= direction[c]
        pts.append(x / (1 << BITS))
    return pts
# real direction numbers: primitive polynomials over GF(2)

python
from scipy.stats import qmc

# 4096 scrambled Sobol points in 20-D
qmc.Sobol(d=20, scramble=True).random_base2(m=12)

Chapter 5: Koksma-Hlawka — error ≤ variation × discrepancy

You finally have a guarantee, not just a hope. For any reasonably smooth integrand, the error of your QMC estimate is at most V(f) times D*_N: the function's "wiggle budget" multiplied by your points' non-uniformity. Halve the discrepancy and you halve the worst-case error — deterministically, no probability hand-waving.

This is the Koksma-Hlawka inequality, the theorem that justifies the whole field. It factors the error cleanly into two independent pieces: a property of the function (its total variation V), and a property of the points (their star discrepancy D*).

The total variation V(f) is the integrand's "wiggle budget" — in 1-D, the sum of all its up-and-down moves, Σ|f(x_i+1) − f(x_i)| in the limit. (In higher dimensions it is the Hardy-Krause variation, a sum of variations of all the mixed partials.) A flat function has V = 0; a wildly oscillating one has large V.

The crucial consequence: once you have fixed the integrand, V(f) is a constant you cannot change. The only lever left is D*_N. That is why QMC research is obsessed with discrepancy — it is the entire controllable factor in the worst-case error. Build lower-discrepancy points and you directly tighten the bound.

Let us see the bound work end to end. Integrate f(x) = x over [0,1] (true value 1/2) using the four van der Corput-style midpoints {1/8, 3/8, 5/8, 7/8}. This f is monotone increasing, so its total variation is simply V(f) = f(1) − f(0) = 1 − 0 = 1. The point set's star discrepancy we computed in Chapter 1: D*_N = 1/(2N) = 1/8 = 0.125.

The QMC estimate is the plain average of the function at those points: (1/8 + 3/8 + 5/8 + 7/8)/4 = (16/8)/4 = 2/4 = 0.5. That equals the true integral exactly, so the actual error here is 0. The Koksma-Hlawka bound is V(f) · D*_N = 1 × 0.125 = 0.125. The bound (0.125) correctly upper-bounds the true error (0) — and it is loose, as promised: the real error sat far below the worst case.

python
def total_variation_1d(f, grid=10000):
    xs = np.linspace(0, 1, grid)
    return np.sum(np.abs(np.diff(f(xs))))   # Hardy-Krause TV

def kh_bound(f, pts):
    V = total_variation_1d(f)
    Dstar = star_discrepancy_1d(np.sort(pts))
    return V * Dstar                       # |error| <= V * D*

kh_bound(lambda x: x, np.array([1,3,5,7])/8)  # 1 * 0.125 = 0.125

python
from scipy.stats import qmc
# no single call gives V(f); pair the
# discrepancy with your own TV estimate
qmc.discrepancy(pts.reshape(-1,1))  # the D* half

Chapter 6: The Rate — (log N)^d/N beats 1/√N

Plain Monte Carlo is brutally honest about cost: to add one decimal digit of accuracy you need 100× more samples, because error falls only as 1/√N. QMC changes the deal. With error decaying nearly as 1/N, one more digit costs only ~10× more points. Over a long simulation that is the difference between an hour and a week.

You already proved the 1/√N law in Monte Carlo — this chapter does not re-derive it, it contrasts against it. The new fact is the QMC rate: a low-discrepancy sequence achieves star discrepancy of order (log N)^d/N, and by Koksma-Hlawka that is the error rate. Strip the slowly-growing (log N)^d factor and the core is 1/N, twice as fast (in the exponent) as 1/√N.

Let us compare the two formulas directly at N = 64, d = 1 — no random draws, just arithmetic. The MC scale is 1/√64 = 1/8 = 0.125. The QMC discrepancy scale with d = 1 is (log N)¹/N = log(64)/64. Now log(64) = 4.1589, so log(64)/64 = 4.1589/64 = 0.06498 — already smaller than the MC scale, even at this modest N.

The deeper point is the cost-to-improve, not the value at one N. To halve MC error from 0.125 to 0.0625, you need 1/√N = 0.0625, i.e. √N = 16, so N = 256 — four times the points. To halve QMC error (ignoring the slow log factor) you roughly double N. So MC's cost-to-improve is quadratically worse: 4× vs 2× for every halving. That asymmetry, compounded over many digits, is the entire reason QMC exists.

But read the formula's fine print. The (log N)^d factor explodes with dimension d. For large d and modest N, (log N)^d can exceed √N, so QMC can actually lose to MC until N is astronomically large. The clean QMC win is for small-to-moderate dimension — which is exactly the cliffhanger the next chapter resolves with "effective dimension."

The canvas draws the convergence on log-log axes, where a power law N^−p appears as a straight line of slope −p. Random MC has slope ≈ −0.5 (and is noisy); Sobol QMC has slope near −1; a regular grid has slope ≈ −1/d (catastrophic past 2-D). Drag the dimension slider and watch the QMC-beats-MC crossover slide rightward as d climbs.

python
def conv_curves(f, true, d, Ns):
    from scipy.stats import qmc
    mc, qm = [], []
    for N in Ns:
        mc.append(abs(np.mean(f(np.random.rand(N, d))) - true))
        X = qmc.Sobol(d=d, scramble=True).random(N)
        qm.append(abs(np.mean(f(X)) - true))
    return np.array(mc), np.array(qm)
# log-log slopes: MC ~ -0.5, QMC ~ -1.0 for small d

python
from scipy.stats import qmc

# average f over these for ~1/N convergence
qmc.Sobol(d=d, scramble=True).random(N)

Chapter 7: Effective Dimension — the rescue from the curse

You face a 360-dimensional integral — the monthly cash flows of a 30-year mortgage pool. The (log N)³⁶⁰ factor in the QMC bound is astronomically large; the theory from last chapter says QMC should be useless. Yet in the famous Paskov-Traub experiment, QMC crushed Monte Carlo on exactly this problem. The resolution is that the integral only looks 360-dimensional.

The escape hatch is effective dimension: most real high-dimensional integrands concentrate almost all their variance in a few "important" coordinates. The nominal dimension d can be 360, but the truncation effective dimension — the smallest k of coordinates that together carry, say, 95% of the variance — may be just 1 or 2.

This rescues QMC because Sobol stratifies its leading coordinates best. The first dimension of a Sobol sequence is the perfectly-balanced van der Corput sequence; later dimensions are progressively less perfectly balanced. So if you arrange for the dominant variance to live in coordinates 1, 2, 3, Sobol equidistributes exactly the directions that matter, and the (log N)³⁶⁰ worst case never materializes.

Let us make this concrete with an integrand dominated by its first coordinate: g(x) = exp(x₁) + 0.05·Σ_i=2..10 sin(x_i) on [0,1]¹⁰. The first term carries almost all the action; the nine small sin terms are perturbations. Its nominal dimension is 10, but its effective dimension is about 1.

Compute the true value by integrating term by term. The first term: ∫₀¹ e^x dx = e − 1 = 1.71828. Each sin term: ∫₀¹ sin(x) dx = 1 − cos(1) = 0.45970. There are nine of them, each scaled by 0.05, contributing 9 × 0.05 × 0.45970 = 0.20687.

Add them: the true value is 1.71828 + 0.20687 = 1.92515. Because the nominal d = 10 but the effective dimension is ~1, Sobol — which stratifies coordinate 1 superbly — integrates g far more accurately than the (log N)¹⁰ worst-case bound would ever suggest. The canvas shows how a Brownian-bridge construction deliberately loads variance into the low coordinates.

python
def brownian_bridge(U, T=1.0):
    from scipy.stats import norm
    Z = norm.ppf(U)            # Gaussian increments via inverse CDF
    # fill endpoint first, then midpoints recursively
    return build_bridge(Z, T)     # leading coords set coarse shape
# effective dim = smallest k explaining ~95% of variance

python
from scipy.stats import qmc, norm

# Gaussian QMC increments for a 360-step path
norm.ppf(qmc.Sobol(d=360, scramble=True).random_base2(m=14))

Chapter 8: Randomized QMC — error bars without losing the rate

Your boss wants a confidence interval. Plain Sobol gives one deterministic number with no error bar — you cannot say "±0.3%". The fix is to randomize the point set in a way that preserves its low discrepancy: apply a random scramble or shift, repeat with a few independent randomizations, and take the sample standard deviation of the resulting estimates as your error bar.

This is randomized QMC (RQMC). The trick is to choose a randomization that keeps the net structure intact. Two standard ones: a Cranley-Patterson shift (add a single random offset to every point, modulo 1) and Owen scrambling (randomly permute the digits of each coordinate). Both leave the points just as evenly spread, so each randomized estimate still converges at the QMC rate — but now it is also unbiased and you can estimate its variance.

Let us trace a Cranley-Patterson shift by hand. Fix the random offset U = 0.8 for this replicate, and apply it to two van der Corput points. Shift 0.375: frac(0.375 + 0.8) = frac(1.175) = 0.175. Shift 0.625: frac(0.625 + 0.8) = frac(1.425) = 0.425.

Notice what the shift did. The two points moved by the same amount and wrapped around the torus (the [0,1] interval glued end to end), so their relative spacing — and hence their low discrepancy — is unchanged, while their absolute positions are now random. Averaging the integrand over many independent shifts is unbiased, and the spread across shifts is your error bar.

Scrambling is unbiased too, and we can confirm it numerically with a deterministic test. Integrate ∏_i(2x_i) over [0,1]² (true value 1.0, since ∫2x dx = 1 in each dimension). Average the estimate over 500 independently-scrambled Sobol sets of 256 points each: the mean comes out to 1.000. No bias — each scramble is centered on the truth, and the scatter of the 500 estimates is the variance you report.

The payoff: a handful of replicates R (say 10–30) buys you an honest standard error std/√R, while each replicate individually still rides the near-1/N QMC convergence. You pay a small multiplier in total cost to gain a confidence interval — far cheaper than plain MC's wide intervals for the same accuracy. The canvas builds the histogram of R estimates live.

python
def rqmc_estimate(f, d, N, R, rng):
    from scipy.stats import qmc
    base = qmc.Sobol(d=d, scramble=False).random(N)
    ests = []
    for _ in range(R):              # R independent randomizations
        U = rng.random(d)
        shifted = (base + U) % 1.0   # Cranley-Patterson shift
        ests.append(np.mean(f(shifted)))
    ests = np.array(ests)
    return ests.mean(), ests.std(ddof=1)/np.sqrt(R)  # est, std err

python
from scipy.stats import qmc

# R Owen-scrambled replicates -> mean +/- std error
[qmc.Sobol(d=d, scramble=True, seed=r).random(N)
 for r in range(R)]

Chapter 9: Showcase — pricing a high-dimensional option

Price an arithmetic Asian option: its payoff depends on the average of a stock's price over many monitoring dates, so each price path is a high-dimensional point. Monte Carlo converges painfully slowly. Watch scrambled Sobol — fed through a Brownian bridge that loads the path's gross shape into its first coordinates — reach the same accuracy with a fraction of the paths, and with honest error bars from randomization.

This is the showcase that made QMC famous. The deterministic anchor is the Paskov-Traub benchmark: they priced a Collateralized Mortgage Obligation with 30 years of monthly cash flows. That is 30 × 12 = 360 monitoring dates, hence a 360-dimensional integral — the dimension that made their result surprising, because the (log N)³⁶⁰ bound predicted disaster.

Why Sobol and not Halton here? Halton's base for the 30th dimension would be the 30th prime, 113, whose radical inverse cycles so slowly it degrades to a near-collinear streak (Chapter 3). Sobol stays in base 2 as a digital net, so it keeps its low discrepancy in all 360 dimensions. And the cost asymmetry to gain one digit — MC needs 100× the paths, QMC only ~10× — compounds over those dimensions into the large practical speedups the dashboard shows.

Recall the halving cost from Chapter 6, now in finance terms: for QMC to halve its error costs about a 2× increase in paths; for MC, about 4× (since (2)² = 4). Over a long pricing run that turns a week of compute into an afternoon. The dashboard below simulates GBM price paths under each method and tracks the running price, its shrinking confidence band, and a "paths to reach target accuracy" counter.

The mechanics in the dashboard: each path is a sequence of log-returns. MC draws them from independent Gaussians; QMC draws Sobol points, maps them to Gaussians via the inverse normal CDF, and (optionally) reorders them through a Brownian bridge so the path's coarse shape rides Sobol's best coordinates. The payoff is max(average price − strike, 0), discounted. Toggle the construction and watch the effective-dimension effect from Chapter 7 sharpen or blunt the QMC advantage.

Play with it: raise the number of monitoring dates (the dimension), lower the path count, switch MC → plain-Sobol → scrambled-Sobol+bridge, and read the speedup ratio (MC paths / QMC paths for equal accuracy). The bridge switch is the most instructive — with it on, QMC's error curve steepens dramatically even at high nominal dimension.

Chapter 10: Showcase — less-grainy rendering with blue-noise samples

Two renders of the same soft shadow, four samples per pixel each. The random-sampled one is splotchy — bright and dark clumps that look like film grain. The low-discrepancy one is smooth, almost noise-free. You spent the same compute. The only difference is where the four samples landed inside each pixel, and that is pure quasi-Monte Carlo.

Every pixel in a ray-traced image is itself an integral: the average of incoming light over the pixel's area, the lens aperture, the light source, and time. A renderer estimates that integral with a few samples per pixel. Random samples clump, so neighboring pixels get wildly different estimates — that variance is the visible grain.

Low-discrepancy samples fix this two ways. First, within a pixel they cover the integration domain evenly, so each pixel's estimate is closer to the truth (lower variance, by exactly the QMC argument). Second, and subtler: good samplers make the residual error blue noise — error energy pushed to high spatial frequencies, which the human eye tolerates far better than low-frequency blotches. The same RMSE looks dramatically cleaner.

Let us anchor the idea with the area integral a single pixel computes. Estimate the area of a quarter-disk inside the unit square (true value π/4 = 0.7854) — a stand-in for a pixel's coverage integral — using the first four 2-D Halton points (bases 2 and 3). A point counts as "inside" if x² + y² ≤ 1.

The first four Halton points are (1/2, 1/3), (1/4, 2/3), (3/4, 1/9), (1/8, 4/9). Check each against the circle: (0.5, 0.333) gives 0.361 ≤ 1, inside; (0.25, 0.667) gives 0.507, inside; (0.75, 0.111) gives 0.575, inside; (0.125, 0.444) gives 0.213, inside. All four are inside, so the estimate is 4/4 = 1.0 against the true 0.7854.

With only four points the estimate is rough — but crucially the Halton points are evenly spread, no two in the same pixel-quadrant. As samples grow, the estimate converges far faster and far smoother than random, with no clumps to create grain. That smoothness, sample after sample, is the anti-grain effect the sandbox below makes visible.

Chapter 11: Connections — where QMC goes next

You have learned to beat the 1/√N wall not by being cleverer about randomness but by abandoning it for structure: deterministic low-discrepancy points that fill space. That same idea now shows up far beyond integration — in how modern probabilistic ML draws its samples — and it stops cleanly where a different idea, reweighting, takes over.

The one sentence that defines QMC: it reduces error by choosing better points on a known domain (usually the unit cube), not by reweighting samples and not by building a chain. Hold that next to the family's other tools and the division of labor is crisp.

python
import torch
from scipy.stats import qmc

a = qmc.Sobol(d=2, scramble=True).random(1024)
b = torch.quasirandom.SobolEngine(
        dimension=2, scramble=True).draw(1024)
# QMC for VI: z = mu + sigma * Phi^-1(sobol)
# better POINTS, not weights -> lower-variance ELBO grads

python
import torch

# scrambled Sobol in the ML stack
torch.quasirandom.SobolEngine(
    dimension=d, scramble=True).draw(N)