An error bar for points that are not random

Two lines of argument about sampling have pulled in opposite directions. The error that does not care how many dimensions established the virtue of random sampling: its error falls like $N^{-1/2}$ whatever the dimension, and the sample itself says how large that error is, because the spread of the sampled values estimates it. Points too even to be random and the lattice rules established the virtue of deterministic points: spread evenly, or placed so their blind spots miss the integrand, they integrate far more accurately — and say nothing at all about how accurately.

That silence is a practical disaster. An estimate without an error bar cannot be trusted, cannot be compared with another, and cannot tell its user when to stop adding points. The deterministic error bounds exist, but they involve quantities like the integrand’s total variation that are harder to compute than the integral itself, and they are usually pessimistic by orders of magnitude.

Randomised quasi-Monte Carlo resolves the dilemma, and it does so with a trick that sounds too simple to work. Take the deterministic point set and randomise it — slightly, carefully, in a way that keeps its evenness — and then repeat the randomisation a few times. Each randomised set gives an unbiased estimate; the repeats are independent; their spread gives a confidence interval exactly as random sampling does. Nothing about this requires the randomisation to be large. The points move, but their evenness does not change: the structure that made them accurate is preserved exactly, and only the part of the construction that was arbitrary anyway — where the whole pattern sits, which way the digits run — is left to chance. And for the best randomisations, the accuracy is not just kept but improved.

Shifting a lattice by a random vector #

The simplest randomisation was proposed by Roy Cranley and Thomas Patterson in 1976. Take the lattice points of a lattice rule, pick one vector $\Delta$ uniformly at random in the unit square, and add it to every point, wrapping round the edges. The shifted set is the same lattice, moved as a rigid whole.

Two facts make this work. First, every shifted point is uniformly distributed over the square — the shift is uniform and wrapping preserves that — so the average of the integrand over the shifted points has expected value exactly equal to the integral. The estimate is unbiased. Second, the shift does nothing to the lattice’s structure: the shifted points are still a coset of the same group, and the Fourier analysis that made the lattice rule accurate applies unchanged. A wave averages to zero over a shifted lattice for exactly the same frequencies as over the unshifted one; only the starting angle of each dual-lattice wave changes. So the shifted rule is exactly as accurate as the original, on average.

Repeat with, say, ten independent shifts. The ten estimates are independent and identically distributed, their mean is an unbiased estimate of the integral, and their sample standard deviation divided by $\sqrt{10}$ is an estimate of its standard error — exactly the statistics of random sampling, applied to ten lattice rules rather than to single points.

The figure uses the same budget for both rows: ten times 233 points. Plain sampling spreads its ten estimates over several hundredths; the shifted lattices, applied to the periodised integrand of the lattice essay, agree to about four decimal places. Both intervals contain the true value, as intervals should about ninety-five times in a hundred. The difference is that one of them is useful.

There is a cost. The ten replicates are ten separate estimates, so the error bar is computed from ten numbers rather than from 2,330, and a confidence interval from ten values is less reliable than one from thousands — the factor $2.26$ in its half-width, rather than $1.96$, comes from treating ten samples honestly. In practice between ten and thirty replicates is the usual compromise: few enough that each replicate has many points and keeps the lattice’s accuracy, many enough that the spread is a trustworthy estimate.

The shift, worked out for one wave #

Why a random shift gives an unbiased estimate, and why its spread measures the lattice’s error, can be seen one wave at a time. Write the integrand as a sum of waves $e^{2\pi i h \cdot x}$, as the lattice essay did. Averaged over the unshifted lattice, a wave gives 1 if its frequency $h$ is on the dual lattice and 0 otherwise. Averaged over the lattice shifted by $\Delta$, the same wave gives $e^{2\pi i h \cdot \Delta}$ if $h$ is on the dual lattice, and still 0 otherwise — the shift multiplies every point’s value of the wave by the same factor.

So the shifted rule’s error is a sum over the nonzero dual frequencies of the integrand’s coefficient times $e^{2\pi i h \cdot \Delta}$. With $\Delta$ uniformly random, each of those factors is a point moving uniformly round the unit circle, with average zero: the error has mean zero, and the estimate is unbiased. And the factors for different frequencies are uncorrelated, so the variance of the error is the sum of the squared sizes of the coefficients on the dual lattice. The variance of a randomly shifted lattice rule is exactly the quantity that measured the unshifted rule’s quality. The replicates’ spread estimates it directly.

That identity is the whole justification of the method. The shift adds no error — the typical size of the shifted rule’s error is the same as the size of the unshifted rule’s — and it converts an unknown, fixed error into a random one whose size can be measured from the data.

Scrambling the digits of a net #

Lattice rules randomise naturally by shifting. The other great family of even point sets, digital nets — of which the Sobol’ sequence is the standard example — are built from binary digits, and they randomise naturally by changing digits.

The two-dimensional Sobol’ sequence has a remarkable property: its first $2^m$ points put exactly one point in every box of area $2^{-m}$ that is formed by halving the square’s sides a whole number of times — every $2^{-a} \times 2^{-b}$ box with $a + b = m$, in every shape from tall and thin to short and wide. That is a much stronger evenness than one point per column and per row.

The simplest randomisation is a random digital shift: pick a random binary string for each coordinate and flip each point’s digits wherever the string has a 1. This is the binary analogue of the lattice’s shift — the points move rigidly in the sense of digit arithmetic — and like the shift it keeps every point uniformly distributed and keeps the box property. Its accuracy, in the opening figure, is about the Sobol’ sequence’s own, falling like $N^{-1}$. That is already a good bargain — the deterministic accuracy, kept whole, with an error bar added — and for lattices the random shift is all that is ever used. For nets it turns out that more randomness buys more accuracy, which is the unexpected part of the story.

Art Owen found a better way in 1995. Nested scrambling flips digits too, but the decision to flip the $k$-th binary digit of a coordinate is made afresh, at random, for each possible value of the first $k - 1$ digits. The first digit is flipped or not for everybody; the second digit is flipped or not separately in the left half and the right half; the third separately in each quarter; and so on down. It is a random permutation of the square that respects its dyadic structure at every scale.

The scrambled points keep the net’s box property exactly — each box of area $2^{-m}$ still holds one point, because scrambling within a box sends it to another box of the same shape — and each individual point is uniformly distributed over the square. What changes is everything the box property does not fix: where inside its box each point sits. A deterministic net places points inside their boxes in a rigid pattern, and the pattern is what limits its accuracy. The scramble randomises the placement independently at every scale.

Why scrambling is better than even spreading #

The surprise in the opening figure is the bottom line: the scrambled points are not merely as good as the deterministic ones, they are better, with error falling like about $N^{-3/2}$ instead of $N^{-1}$.

Owen proved in 1997 that for smooth integrands, the variance of a scrambled-net estimate falls like $N^{-3}$ up to logarithmic factors, so the root-mean-square error falls like $N^{-3/2}$. The reason is cancellation. In a deterministic net, the error comes from the systematic placement of points within their boxes: each box contributes an error of about the same sign, and the errors add. In a scrambled net, the placement within each box is independent and random, so each box’s error has mean zero, and the errors from different boxes partly cancel, like the steps of a random walk. The box structure removes the large-scale error; the randomness turns the small-scale error into noise that averages out.

The digital shift does not do this, because it moves every box’s point in the same way — one random string applied everywhere — so the within-box errors stay correlated from box to box. Nested scrambling makes them independent, and independence is what buys the extra half power.

The picture is worth making concrete. Split the square into $N$ boxes, one point in each. For a smooth integrand, the error a single box contributes depends on where in the box its point sits, relative to the box’s centre of mass for the integrand, and is of order $N^{-1}$ times the box’s area — tiny, but of the same sign in neighbouring boxes if the points sit in the same relative position. Added over $N$ boxes, errors of one sign give a total of order $N^{-1}$, which is the deterministic net’s rate. Make each box’s placement independent and random and the $N$ errors have mean zero; their sum is of order $\sqrt N$ times each one, which is $N^{-3/2}$. That is exactly the gain the opening figure measures, and it is the same arithmetic that makes a random walk wander only as far as the square root of its length.

So the history of these methods can be read as a sequence of exponents. Independent points: $-\tfrac12$. Evenly spread points: about $-1$, and no error bar. Evenly spread points randomised to keep an error bar: about $-1$, with one. Evenly spread points randomised at every scale: about $-\tfrac32$, with one. Each step buys either accuracy or honesty, and the last buys both.

How to spend a fixed budget #

Randomising a point set raises a practical question that plain sampling never faces. With a budget of, say, ten thousand evaluations, should they be spent as ten replicates of a thousand points, or a hundred replicates of a hundred?

For plain Monte Carlo it makes no difference: the error of the combined estimate depends only on the total number of points. For randomised quasi-Monte Carlo it matters a great deal. The error of a single replicate of $N$ points falls like $N^{-\alpha}$ with $\alpha$ between 1 and $\tfrac32$; averaging $R$ independent replicates divides that by $\sqrt R$. With the total $RN$ fixed, the error behaves like $N^{-\alpha + 1/2}$ times a constant — so, whenever $\alpha > \tfrac12$, it pays to make $N$ as large as possible and $R$ as small as the confidence interval allows. That is why the usual practice is between ten and thirty replicates: the fewest that give a trustworthy estimate of the spread, with everything else put into the size of each replicate.

The same reasoning explains why randomised quasi-Monte Carlo combines well with the other variance-reduction tools. Sampling where the answer lives changes the integrand so that it is flatter; flatter integrands have smaller Fourier coefficients and smaller variation, which is exactly what scrambled nets and shifted lattices exploit. The two improvements multiply rather than compete. The one method that does not combine easily is a walk that samples a distribution, whose points are correlated by construction, although versions of quasi-Monte Carlo for Markov chains have been developed.

Where it is used, and why the dimension matters again #

Randomised quasi-Monte Carlo became standard in quantitative finance after 1995, when Spassimir Paskov and Joseph Traub reported that Sobol’ points priced a mortgage-backed security — an integral in 360 dimensions, one per month of a thirty-year mortgage — far more accurately than random sampling. That result was a shock, because the theory said evenly spread points should lose their advantage long before 360 dimensions.

The explanation is effective dimension. The mortgage integrand, like many financial ones, depends mostly on a few coordinates or a few combinations of them, and on the rest only weakly; its effective dimension is small. Quasi-Monte Carlo methods exploit exactly that, and randomisation lets the user see whether the exploitation is working, since the error bar from the replicates reports the actual accuracy rather than a worst-case bound.

That is the deepest benefit of the randomisation. A deterministic quasi-Monte Carlo result in 360 dimensions comes with a bound too pessimistic to use and no other information; a randomised one comes with an interval, and if the effective dimension is too high for the method to help, the interval is wide and says so.

What the pictures cannot show #

The error curves are for one smooth integrand in two dimensions, and the rates they show are properties of that integrand. Owen’s $N^{-3/2}$ needs smoothness; for an integrand with a jump inside the square, scrambled nets improve on plain sampling by a smaller power, and for very rough integrands by little at all. The slopes in the captions are measurements on this function, not guarantees. A reader who wants to know how the methods behave on an integrand of their own has, in the randomised methods, the means to find out: run the replicates and look at the spread. That is a capability the deterministic methods never had, and it matters more in practice than any rate.

The figures also show root-mean-square errors estimated from 24 replicates, which are themselves random; a different set of randomisations would move each point on the error curves a little. And the confidence intervals in the replicate figures rely on the estimates being roughly normally distributed, which is true for large enough point sets and can fail for small ones, where a single replicate can land in an unlucky arrangement.

Still open: the right randomisation, and honest intervals #

The theory of scrambled nets is well developed for smooth integrands, and the practical methods are mature. What remains open is at the edges. For integrands of limited smoothness, the exact rates of scrambled nets and shifted lattices are known only in special classes, and which randomisation is best for a given kind of irregularity is not settled.

The reliability of the intervals is a second open question. Because the variance of a scrambled estimate falls so fast, the distribution of the estimate is not always close to normal at practical sizes, and intervals built from a small number of replicates can be too narrow. How many replicates are needed for an interval to have its stated coverage — and how to build intervals that are honest when the estimates are skewed — is an active topic, and an important one, because an error bar that is too small is worse than none. Bootstrap intervals, intervals built from the order statistics of the replicates, and intervals that use a known bound on the estimate’s skewness have all been proposed, and each trades width for reliability in a different way; none is established as the right default.