Points too even to be random

Independent random points look, to the eye, badly distributed. They arrive in clumps with gaps between, and a scatter of two hundred of them across a square will have a knot in one corner and an empty patch in another. That is not a failure of the generator; it is what independence produces, and anybody laying out points by hand would do considerably better.

The error of a Monte Carlo estimate falls like $1/\sqrt N$ precisely because of that clumping. The question this rung asks is what happens if the clumping is simply prevented.

Measuring evenness #

Before anything can be improved, evenness has to be a number. The standard one is discrepancy: over all boxes with one corner at the origin, the largest difference between the share of the points inside and the area of the box.

$$D_N = \sup_B \left|\frac{\#\{\text{points in }B\}}{N} - \operatorname{area}(B)\right|.$$

A perfectly even set has small discrepancy; a set with a clump has a box around the clump holding more than its share, so its discrepancy is large.

The figures measure this over a grid of test boxes rather than over all of them, which gives a lower bound on the true supremum — the worst box has a corner at a data point and the grid need not contain it. That is stated in the caption, because it matters: a lower bound is enough to establish that one set is worse than another, which is what the comparison needs, and is not enough to certify that any set is good.

There is a second reason to define evenness by boxes with a corner at the origin rather than by all boxes, or by discs, or by anything else. Every such definition gives a slightly different number, and they are all equivalent up to constants in a fixed dimension — so the rate is the same whichever is used, and the rate is what everything here turns on. The corner-anchored version is chosen because it is the one that appears in the error bound below, and because it is the cheapest to compute.

That is worth noticing as a general habit: when several reasonable definitions of a quantity exist and all give the same asymptotics, the one to standardise on is the one a theorem uses.

The rate #

For independent random points the discrepancy is about $\sqrt{\log\log N / N}$, so it falls essentially like $1/\sqrt N$ — the same rate as the integration error, for the same reason.

For a well-chosen deterministic sequence it falls like $(\log N)^d/N$, which for any fixed dimension is very nearly $1/N$. Squaring the rate: a hundred points of such a sequence are worth ten thousand random ones, in a sense the bound below makes precise.

The figure measures both slopes rather than quoting them. It computes the discrepancy of each kind of point set at six sizes from 64 to 2,048, fits a line to the logarithms, and asserts that the random slope is $-0.5$ to within 0.12 and that the Halton slope is steeper by at least 0.2.

How the points are made #

The Halton sequence is built out of the digits of the counting numbers, reflected about the point.

Write $n$ in base $b$. Reverse the digits and put them after a decimal point in the same base. That is the $n$-th term of the van der Corput sequence in base $b$: in base 2 it gives $1/2$, $1/4$, $3/4$, $1/8$, $5/8$, $3/8$, $7/8$, and so on — each new term landing in the largest remaining gap.

That is not a metaphor. The sequence halves the largest gap at every step — the same behaviour that three gaps and no more describes for a different construction, then halves the next largest, and so on, which is exactly the behaviour that keeps the discrepancy small. A $d$-dimensional Halton sequence uses a different prime base for each coordinate, so that the coordinates do not fall into step with one another.

Why reversing the digits works #

The construction looks arbitrary and is not, and the reason is worth a paragraph.

Consider the first $2^k$ terms of the base-2 sequence. Their binary expansions are the reversals of $1$ through $2^k$, which are exactly all the $k$-bit strings — so those $2^k$ points are precisely the points $j/2^k$, one in each interval of length $2^{-k}$. Every dyadic interval of that size contains exactly one point.

That property holds at every scale simultaneously, and it is what small discrepancy means: a box of any size contains its fair share. The reversal is doing something specific — it converts the counting numbers’ habit of filling the low bits first into the habit of filling the finest scale last, so the sequence refines uniformly rather than sweeping from one end.

The same construction in base 3 gives one point in each third, then each ninth, and so on. Using bases 2 and 3 for two coordinates makes them refine at incommensurable rates, so no correlation builds up between them — and using two bases that shared a factor would give points on a lattice rather than a filled square, which is why the bases are primes.

The bound that connects evenness to error #

The link between discrepancy and integration error is the Koksma–Hlawka inequality:

$$\left|\frac1N\sum_{i} f(x_i) - \int f\right| \leq V(f)\, D_N,$$

where $V(f)$ measures how much $f$ wiggles — its total variation in the sense of Hardy and Krause.

Two things about this bound are worth separating.

It is deterministic. There is no probability in it. A Monte Carlo error bar is a statement about a distribution of possible errors; this is a guarantee about the error of the run actually performed. That is a genuinely different kind of statement and is the strongest argument for these methods.

It is nearly useless as a number. Computing $V(f)$ is harder than computing the integral, and computing $D_N$ exactly is itself an expensive combinatorial problem. So the bound is used to establish the rate and essentially never to produce an actual error estimate.

That leaves a practical awkwardness. Quasi-Monte Carlo converges faster and reports no error, while Monte Carlo converges more slowly and reports one. The usual response — randomised quasi-Monte Carlo, in which a low-discrepancy set is randomly shifted and several such shifts are averaged — recovers an error estimate by making the whole thing random again at a coarser level, and keeps most of the rate.

Where the advantage goes #

The $(\log N)^d$ in the numerator is doing real damage in high dimensions and it is honest to say so.

At $d = 30$ the bound $(\log N)^{30}/N$ does not beat $1/\sqrt N$ until $N$ is astronomically large, so the theory says quasi-random points should be useless for the very problems Monte Carlo was invented for. And yet they are used, successfully, on financial problems in hundreds of dimensions.

The resolution is that such problems have a low effective dimension: most of the variance comes from a handful of coordinates, and the low-discrepancy construction is good in exactly those first coordinates. Halton’s later bases are large primes and its higher coordinates are visibly poor, so a construction that puts the important variables first works and one that does not is worse than random. Sobol sequences are used in practice partly because their coordinates degrade far more slowly.

This is one of the places where the theory is pessimistic and the practice is fine, and the honest summary is that the bound is worst-case, most problems are not the worst case, and nobody has a satisfying general account of which are.

What was given up #

The word “random” has left, and several things went with it.

Independence. No error estimate comes free from the run, as noted above.

Robustness to the integrand. A random sample is good against every integrand at once, , because the argument uses nothing about $f$ except its variance. A fixed point set is good against integrands that do not conspire with it — and one can be built that does. Deterministic points have an adversary; random ones do not. It is the same distinction as between a rule that can be gamed and one that cannot, transplanted from voting to numerical integration.

The ability to stop anywhere. Many low-discrepancy constructions are only good at particular sample sizes — powers of two for a Sobol sequence — so stopping after 1,500 points may be worse than stopping after 1,024.

What is gained is the rate and the deterministic bound. Whether the trade is worth making is a question about the problem, and the fact that both methods remain in wide use is the honest indication that neither dominates.

The one-dimensional case, which is a different subject #

In one dimension all of this is older and sharper, and it is worth naming because it connects to a part of this collection that looks unrelated.

The points $\{n\alpha\}$ — the fractional parts of multiples of an irrational $\alpha$ — are equidistributed for every irrational, and how evenly depends on how badly $\alpha$ is approximated by fractions. A number with a very good rational approximation $p/q$ has its first $q$ multiples nearly repeating, which is a clump; a number with no good approximations never does. The golden ratio, the hardest number to approximate, gives the most even sequence, and the gaps between consecutive points take only ever three distinct values.

So the one-dimensional theory of low-discrepancy points is a chapter of the theory of continued fractions, and the quantity governing it is how well a number can be approximated. That is a satisfying place for the subject to come from and it does not generalise: in higher dimensions there is no analogue of continued fractions with the same power, which is part of why the constructions become combinatorial rather than arithmetic.

What the eye was right about #

It is worth closing the loop on the observation the essay opened with, because it is unusual for an aesthetic judgement to be correct.

A scatter of independent points does look badly distributed, and it is badly distributed by the measure that matters for integration. The clumps are real, the gaps are real, and the reason a hand-placed set looks better is that a person placing points by hand is minimising something like discrepancy without knowing the word.

The usual lesson drawn from random scatter — that people are bad at recognising randomness, and expect it to look more even than it is — is true and is a different point. Both are true at once: a random scatter looks less even than people expect, and is less even than is useful. The first is a fact about intuition and the second is a fact about integration, and confusing them leads to the conclusion that the clumping must be fine because randomness is supposed to look like that.

For deciding whether a pattern is random, the clumping is exactly what to expect. For estimating an integral, it is what to eliminate.

One more comparison is worth making explicit, since it is the reason this rung follows the Metropolis one rather than preceding it. A chain gives up independence in order to sample a distribution that cannot be drawn from; a quasi-random set gives up randomness in order to sample the uniform distribution better than randomness can. Both trades abandon the same property and buy different things, and neither can be applied to the other’s problem — there is no low-discrepancy sequence for an unnormalised target, because evenness is defined against a measure that has to be known.

What the picture cannot show #

The discrepancy printed under each panel is a lower bound, computed over a grid of test boxes. The true supremum is over all boxes and is attained at one whose corner sits at a data point; computing it exactly is possible and expensive, and in two dimensions with a few hundred points the grid is close. In higher dimensions the exact computation is out of reach entirely, so discrepancy — the quantity the whole subject is about — is essentially never measured on the sets that matter.

The pictures are also two-dimensional, and evenness in two dimensions is something the eye judges well. That is why the figures are persuasive and why they are misleading: in twenty dimensions no set of points looks like anything, every pair is nearly the same distance apart, and intuitions formed on a square do not survive — which is the same warning the dimension-counting rung gives about grids, arriving here about the eye instead.

And the $(\log N)^d$ cannot be drawn. It is invisible at the sizes shown, and it is the term that decides whether any of this helps at the dimensions where the methods are used.

The ladder from here #

Rungs above: Sobol and Niederreiter sequences, and the digital-net construction they come from. Randomised quasi-Monte Carlo, and the scrambling that restores an error estimate. Effective dimension and the ANOVA decomposition that measures it. Lattice rules, which are the higher-dimensional relatives of $\{n\alpha\}$. And the lower bounds on discrepancy, which say how even a set of $N$ points can possibly be — a question with a known answer in one dimension and open in general.

Where the two methods meet #

The neatest way to hold both rungs at once is to notice that they are attacking different factors of the same expression.

Monte Carlo error is roughly variability of the integrand times how unevenly the points fall, and the two methods split along that product. Importance sampling attacks the first factor: it does not change how the points clump, it changes the function being averaged so that clumping matters less. Quasi-random points attack the second: they leave the integrand alone and make the points fall evenly.

The two are compatible, and combining them is standard practice — a low-discrepancy set pushed through the inverse of a well-chosen proposal gives points that are both evenly spread and concentrated where the integrand is. Neither idea interferes with the other, because they are about different halves of the same product.

What neither can do is change the exponent from the underlying source: Monte Carlo’s half is a fact about independent averaging, and quasi-random points’ near-one is a fact about how evenly $N$ points can be placed in a cube. Those are the two floors, and everything in this ladder is a matter of getting near one of them.

The trade in one sentence #

Randomness was never the point.

Monte Carlo is called a random method and the randomness is a means: what the argument actually needs is that the sample represents the domain, and independence is one way of arranging that — a cheap, robust, universally applicable way, whose price is the clumping. Choosing points to be representative directly is another way, more effective and less robust.

Naming what a technique actually requires, as against what it happens to use, is worth doing whenever a method is being improved. Here the requirement was evenness, the tool was independence, and the improvement came from noticing that they are not the same thing.