Points on a lattice that see almost nothing

Points too even to be random showed that replacing random sample points by a deterministic sequence spread as evenly as possible — Halton’s, built from digits in two bases — improves the error of numerical integration from about $N^{-1/2}$ to about $N^{-1}$, at the cost of the error bar that random sampling provides for free. It ended by naming the constructions that go further, and one of them does something quite different from spreading points evenly.

A lattice rule puts its points on a lattice, tilted and wrapped round the unit square. The points are not especially even by the measures that judge Halton’s sequence. But for a large class of integrands they are extraordinarily accurate — their error can fall like $N^{-2}$, or faster — and the reason is not geometry but Fourier analysis. A lattice rule sees almost every frequency perfectly, and is completely blind to the rest. The idea goes back to Nikolai Korobov and Edmund Hlawka around 1960, and it is a clean example of a method whose power comes from knowing exactly where it fails: the failures are listed in advance, as a lattice of frequencies, and the design problem is to push that list out to where the integrand has nothing to lose.

A lattice wrapped round the square #

A rank-one lattice rule with $N$ points and generator $(1, z)$ uses the points

$$\left(\frac{k}{N}, \left\{\frac{kz}{N}\right\}\right), \qquad k = 0, 1, \ldots, N - 1,$$

where $\{\cdot\}$ takes the fractional part: walk in a straight line of slope $z$ and wrap round whenever the walk leaves the square. With $N$ and $z$ chosen well the points fill the square in a regular tilted grid. With them chosen badly — $z = 1$, say — every point lands on the diagonal, and the rule is useless; the whole art is in the choice of $z$.

The best-known choice in two dimensions takes $N$ to be a Fibonacci number and $z$ the one before it — $N = 89$, $z = 55$, for instance. These are the Fibonacci lattices, and they are the best two-dimensional lattice rules in a sense made precise below. Their connection to Fibonacci numbers is the same one the fractions that beat every smaller one found: $55/89$ is a continued-fraction convergent of the golden ratio’s reciprocal, the number worst approximated by fractions, and that is exactly what keeps the lattice’s points from lining up on too few lines.

The lattice has an obvious evenness — one point in every column and row — and an obvious regularity: it lies on families of parallel lines. That regularity is the thing a recurrence that stands in for chance was ashamed of, since points on few lines make bad random numbers. For integration it is an asset, and the reason needs a change of viewpoint.

It helps to see the lattice as a group. Adding two lattice points and wrapping round gives another lattice point: the $j$-th plus the $k$-th is the $(j + k)$-th. So the point set is closed under addition modulo 1, a finite subgroup of the torus, and averaging over it is averaging over a group. Everything that follows — the exactness, the blind spots — is the standard behaviour of a character averaged over a finite group: it averages to zero unless it is trivial on the group, and then it averages to one. Characters averaged over residue classes did the same job for primes in arithmetic progressions; here the group is a set of points in a square and the characters are waves.

What a lattice rule cannot see #

Write the integrand as a Fourier series — a sum of waves $e^{2\pi i (h_1 x + h_2 y)}$ over integer frequencies $(h_1, h_2)$, as a square wave is built from round ones. The integral over the square is the coefficient of the constant wave, $(0, 0)$; every other wave integrates to zero. So a rule is exact for a function if it averages every non-constant wave to zero.

Now average a single wave over the lattice points. The wave at the $k$-th point is $e^{2\pi i k (h_1 + h_2 z)/N}$, a power of one complex number, and the average of $N$ consecutive powers of a root of unity is zero — unless the root is $1$, when the average is $1$. That happens exactly when $h_1 + h_2 z$ is divisible by $N$.

So the rule averages every wave to zero except those whose frequency satisfies $h_1 + z h_2 \equiv 0 \pmod N$. Those frequencies form a lattice of their own, the dual lattice, and on them the rule reads the wave as a constant.

The error of the lattice rule is therefore exactly the sum of the integrand’s Fourier coefficients at the nonzero points of the dual lattice. Nothing else contributes. That changes the design problem: to make the rule accurate for smooth integrands, whose Fourier coefficients fall off rapidly as frequencies grow, choose the lattice so that its dual has no small nonzero points.

“Small” here has a particular meaning. For a smooth function of two variables, the coefficient at $(h_1, h_2)$ is roughly of size $1/(\bar h_1 \bar h_2)^\alpha$, where $\bar h = \max(1, |h|)$ and $\alpha$ measures smoothness. So the dual points that matter are the ones with small $\bar h_1 \bar h_2$, and the quality of a lattice is the smallest such product over its dual. For Fibonacci lattices that product grows in proportion to $N$, the best possible, and it is why the error falls like $N^{-\alpha}$ up to a logarithm: for the integrand in the opening figure, $\alpha = 2$.

One wave, averaged by hand #

The mechanism is simple enough to check without the figure. Take the 34-point Fibonacci rule, with generator $z = 21$, and the wave with frequency $(1, 1)$, which oscillates once across the square in each direction. At the $k$-th lattice point the wave has turned through $k(1 + 21)/34 = 22k/34$ of a full turn. As $k$ runs from 0 to 33 that angle takes 34 values spaced evenly round the circle — since 22 and 34 share only the factor 2, it goes round twice through 17 values each — and the average of the wave over them is exactly zero. The rule gets the integral of this wave right: zero.

Now take the frequency $(-13, -1)$, the ringed point nearest the origin in the dual-lattice figure. The angle turned at the $k$-th point is $k(-13 - 21)/34 = -k$ full turns: a whole number, so the wave equals 1 at every lattice point, and the rule reports its average as 1 when the true integral is 0. Every frequency either behaves like $(1, 1)$ — averaged to zero, perfectly — or like $(-13, -1)$ — mistaken for a constant, completely. There is nothing in between. The rule’s accuracy on a given integrand is therefore exactly the question of how much of the integrand lives at the second kind of frequency.

For the smooth integrand in the opening figure, the coefficient at $(-13, -1)$ is about $1/(13^2 \cdot 1^2)$ times a constant, a little under one per cent of the leading term, and the coefficients at the other dual points are smaller still. That is the whole error. Double the number of points to the next Fibonacci lattices and the nearest dual point moves out to products around 34, then 89, and the error falls as the square of that distance.

Why the rate beats even spreading #

The Halton sequence is judged by discrepancy — how evenly it fills every box — and the bound that turns discrepancy into error gives about $N^{-1}$ for any integrand of bounded variation. That bound is uniform over a huge class of functions, and it is sharp for that class: some function of bounded variation really does have error of order $1/N$ against a set of low discrepancy.

A lattice rule plays a different game. It gives up uniform performance and targets smooth periodic functions, for which the Fourier coefficients decay fast, and it arranges for its blind spots to lie where those coefficients are negligible. The payoff grows with the smoothness: twice-differentiable periodic integrands get about $N^{-2}$, smoother ones faster still. In the opening figure the lattice’s slope is about $-1.8$ against Halton’s $-1$, and by 1,597 points the lattice is two orders of magnitude ahead of Halton and four ahead of random sampling. The comparison is not a contest between two answers to one question. Halton’s sequence answers “how evenly can points be spread?” and a lattice rule answers “where should the unavoidable blind spots go?” — and for smooth periodic integrands the second question has the better answer.

When the integrand does not wrap round #

The analysis above assumed a periodic integrand — one whose values and slopes on the right edge of the square match those on the left, and top matches bottom. Most integrands in practice are not like that, and the difference is decisive.

A lattice rule sees the square as a torus, with opposite edges glued. An integrand that does not match across the edges has, on that torus, a discontinuity along the seam — and a discontinuous function has Fourier coefficients that decay only like $1/h$, slowly enough that the dual lattice’s points pick up substantial error. The lattice’s advantage collapses to roughly the Halton rate, as the figure shows.

The repair is a change of variable. Substitute $x = \phi(t)$ in each coordinate, with $\phi$ an increasing map of $[0, 1]$ onto itself whose derivative vanishes at both ends — the figure uses $\phi(t) = 3t^2 - 2t^3$, with derivative $6t(1 - t)$. The integral is unchanged, since it is just rewritten in new coordinates, but the new integrand, multiplied by the derivative, vanishes at the edges together with its slope. It is effectively periodic.

Periodising is the standard way lattice rules are applied to real integrands. It costs one change of variable, and the lattice’s rate returns. It does nothing for random points, whose error depends only on the variance of the integrand, or for the Halton sequence, whose error depends on its variation — the change of variable helps only a method whose blind spots are Fourier frequencies.

A generator that is also a number-theory problem #

Choosing a good lattice is a question about the generator $z$, and it is a number-theory question. The dual lattice has a short vector exactly when $z/N$ is well approximated by fractions with small denominators — when $h_1 + z h_2 \equiv 0 \pmod N$ has a solution with small $h_1, h_2$. So the best generators are those for which $z/N$ is badly approximable: its continued fraction has small partial quotients.

That is why Fibonacci ratios are best in two dimensions, since their continued fractions are all ones, and it is the same question Stanisław Zaremba asked in 1971 — whether every denominator $N$ has some $z$ with all partial quotients of $z/N$ at most 5 — which the essay on the Stern–Brocot tree’s fans recorded as open. Zaremba asked it because of lattice rules. In higher dimensions there is no continued fraction, and good generators are found by search; Ian Sloan and his collaborators showed in the early 2000s that choosing the generator’s components one at a time, each to minimise a computable error bound given the ones before, produces lattices within a constant of the best possible — the component-by-component construction now used in practice. The construction is greedy in the same sense as greedy colouring: it never revisits an earlier choice. Greedy colouring can be arbitrarily bad, depending on the order; the surprise is that for lattice rules, greedy is provably almost as good as optimal, because the error bound being minimised splits into a sum over coordinates in which each new component’s contribution can be controlled given the old ones. The search for each component is over at most $N$ candidates, so a lattice in a thousand dimensions with a million points is found in minutes.

Where lattices earn their keep #

In a handful of dimensions a lattice rule competes with methods that are easier to use. Its real territory is high dimension, where the alternatives collapse. A product of one-dimensional rules with ten points in each of fifty coordinates needs $10^{50}$ evaluations; a lattice rule with a million points is a million evaluations whatever the dimension. The error that does not care how many dimensions is random sampling’s great virtue, and a lattice rule can keep that indifference while improving the rate — provided the integrand’s dependence on its many coordinates is uneven.

That proviso was made precise by Ian Sloan and Henryk Woźniakowski in 1998. They weighted the coordinates by importance and showed that lattice rules achieve errors independent of the dimension when the weights shrink fast enough — when the integrand depends strongly on a few coordinates and weakly on the rest, which is what many integrands in finance and physics do. The component-by-component construction then finds a generator tuned to the weights.

The largest recent application is in uncertainty quantification: computing the expected behaviour of a physical system — groundwater flowing through rock of uncertain permeability, say — whose inputs are random fields described by hundreds or thousands of random parameters. Each evaluation of the integrand is the solution of a differential equation, far too expensive to repeat millions of times, and lattice rules with weights derived from the equation have been shown, by Frances Kuo, Christoph Schwab and Sloan in 2012, to reach a given accuracy with far fewer solves than random sampling needs.

What the lattices cannot show #

The figures are two-dimensional, which is where lattice rules are least needed: in two dimensions tensor-product rules and adaptive methods are hard to beat. The reason lattice rules matter is high dimension, where their error bounds, suitably weighted, avoid the exponential dependence on dimension that grids suffer — and nothing drawn in a square shows that.

The error curves are also for integrands chosen to be smooth, and the rates in the captions are properties of those integrands. An integrand with a kink or a jump inside the square defeats every lattice rule, periodised or not, since its Fourier coefficients decay slowly in directions no change of variable can repair. And the curves say nothing about the error on a particular problem in advance, because the lattice rule, being deterministic, gives no estimate of its own error: that is the price the previous comparison of even and random points identified, and it is the problem the essay on randomising these points takes up. A reader looking at the lattice’s steep line should also remember that its steepness is a statement about many values of $N$: at any single $N$ the lattice gives one number and no indication of how far it is from the truth.

Still open: the best lattices in high dimension #

In two dimensions the best lattice rules are known — the Fibonacci lattices — and their error is understood up to constants. In higher dimensions the picture is incomplete. The component-by-component construction gives lattices that achieve the optimal rate in the worst case over a weighted class of functions, but “optimal” there means up to constants and logarithmic factors, and the constants can grow with the dimension in ways that are known only for particular weightings.

The underlying number-theory question is open in every dimension above two: for each $N$, what is the lattice whose dual has the largest minimum product $\bar h_1 \bar h_2 \cdots \bar h_s$, and how does that minimum grow with $N$? The answer in two dimensions is Zaremba’s question in another form, and it is not settled; in higher dimensions even the order of growth of the best possible minimum is known only within logarithmic factors. And whether lattices are the right structure at all is not settled either: digital nets, built from binary digits in the way the Sobol’ sequence is, compete with lattices on smooth integrands, and for some weighted classes of functions it is not known which achieves the smaller constant.