The error that does not care how many dimensions

Dropping needles to estimate pi is a bad way to compute pi and a good way to introduce a method, and the essay that introduced it said so: the integral being estimated is one a row of rectangles would do far better on, and far faster.

That is true in one dimension and it stops being true, sharply, somewhere around four. This rung is about where the crossover is and why it is there.

Two ways to average a function #

An integral over a cube is an average — the average value of the function, times the volume. Two ways of estimating an average suggest themselves.

Take a grid. Divide each axis into $m$ pieces, evaluate the function at the middle of each cell, and average. That uses $N = m^d$ points in $d$ dimensions. For a smooth function the midpoint rule’s error in one dimension falls like $m^{-2}$, and the product rule inherits it, so the error is $N^{-2/d}$.

Take random points. Pick $N$ points uniformly in the cube, evaluate, and average. The law of large numbers says this converges; the central limit theorem says the error is about $\sigma/\sqrt N$, where $\sigma$ is the standard deviation of the function’s values. There is no $d$ anywhere in that expression.

Those two rates are the whole subject. One improves fast and degrades with dimension; the other improves slowly and does not degrade at all.

Where they cross #

Setting $N^{-2/d}$ against $N^{-1/2}$ gives the crossover at $d = 4$. Below four dimensions the grid wins and keeps winning; above four the random points win, and the margin grows without limit.

In eight dimensions the grid’s error falls like $N^{-1/4}$, so improving it by a factor of ten costs ten thousand times as many points. The random average improves by a factor of ten for a hundred times as many. In twenty dimensions the grid needs $10^{20}$ points for one significant figure and the random average needs a hundred for the same.

The figure computes each grid’s slope by least squares over its own points, and asserts that the measured slope is $-2/d$ to within 0.06. The random slope is measured the same way from the spread of twenty-four independent runs at each sample size, and asserted to be $-0.5$ to within 0.09. Neither number is drawn from the formula it is being compared with.

How a grid of a hundred million points is summed #

There is a difficulty in drawing the figure at all that is worth admitting, because the way round it is also an argument.

A grid in eight dimensions with ten points per axis has $10^8$ cells, which is a hundred million function evaluations for a single point on a single line. Several such points, for several dimensions, is not a figure anybody can generate on a page.

The integrand is therefore chosen to be a product — the same bump on every axis — so that the integral over the cube is the one-dimensional integral raised to the power $d$, and the product grid’s estimate is likewise the one-dimensional rule raised to the power $d$. That is not an approximation to the grid; it is exactly what the grid gives, and it takes $m$ evaluations rather than $m^d$.

The caption says so, because a reader is entitled to ask how $10^8$ points were summed and the answer they were not, and did not need to be is part of the content. It also shows what makes high-dimensional integration hard: the difficulty is not the arithmetic of any one evaluation, it is that a grid demands exponentially many of them, and a product integrand is exactly the special case where that demand can be evaded.

For a genuinely non-separable integrand there is no such evasion, and the exponential cost is real. The random estimate does not care either way.

What the crossover really depends on #

The number four is not a constant of nature and it is worth taking apart, because two of the three ingredients could have been different.

The grid’s exponent $2/d$ comes from the midpoint rule’s order of accuracy, which is two. A better one-dimensional rule changes it: Simpson’s rule is fourth order, giving $N^{-4/d}$ and a crossover at eight; a Gaussian rule of high order does better still, at least for very smooth integrands. So four is the crossover for the crudest reasonable grid, and better grids push it up.

The random exponent is a half and nothing improves it, so the crossover moves in one direction only.

What does not move at all is that a crossover exists. Any product rule has an exponent of the form $c/d$ for some fixed $c$, and $c/d$ falls below $1/2$ once $d$ exceeds $2c$. There is no quadrature rule of any order for which the grid keeps winning as the dimension grows, because the exponent is always divided by $d$ and the random one never is.

That is the content of the phrase curse of dimensionality, and it is worth having in this precise form. It does not say high dimensions are hard; it says any method whose cost is a product over the axes has an exponent that dies like $1/d$, and any method that does not build a product does not.

The variance is where the dimension hides #

The claim that Monte Carlo error is dimension-free needs one qualification, and it is the qualification that decides whether the method is useful in practice.

The error is $\sigma/\sqrt N$, and $\sigma$ is the standard deviation of the function over the cube. That constant can depend on dimension, and often badly. A function that is a product of $d$ bumps has a standard deviation that grows or shrinks exponentially in $d$ relative to its mean, so the relative error at fixed $N$ can be terrible even though the rate is unchanged.

So the honest statement is: the exponent does not depend on the dimension, and the constant does. Everything the method is good for follows from the first half, and everything that makes it hard in practice is in the second — which is why the next rung is about reducing the constant.

What a grid is spending its points on #

There is a way of seeing why the grid fails that has nothing to do with exponents, and it is worth having because it explains the phenomenon rather than measuring it.

A grid with $m$ points per axis places $m^d$ points, and almost all of them are near the boundary of the cube. In eight dimensions with ten points per axis, the fraction of cells with every coordinate in the middle eighth of its range is $(1/8)^8$, about one in sixteen million. The grid has spent essentially its whole budget on corners.

Random points do the same thing — a uniform point in a high-dimensional cube is also near the boundary — but they do not do it systematically. The grid’s failure is that its resolution along each axis is $m = N^{1/d}$, which for $N = 10^6$ and $d = 8$ is about five. Five points per axis is not a resolution; it is barely a sample, and no amount of total points fixes it because the total is being divided among the axes multiplicatively.

A grid resolves each axis to $N^{1/d}$ and a random sample resolves nothing and estimates everything. That is the trade, stated without any exponents, and it is why the crossover is where it is: at $d = 4$, a million points give a grid thirty-one per axis, which is about as coarse as a grid can be and still be worth calling one.

Why the square root, and why it cannot be beaten #

The rate $1/\sqrt N$ is not a property of the algorithm; it is a property of averaging independent things, and it is unimprovable for any method that samples independently.

The variance of an average of $N$ independent draws is the variance of one, divided by $N$; standard deviation is the square root of variance; hence $1/\sqrt N$. Every step is forced, and the same arithmetic governs the wobble around the mean and the width of a random walk after $N$ steps — a walk’s distance from the origin grows like $\sqrt N$ for the identical reason, with the walk’s position playing the part of the accumulated error.

The unimprovability is worth stating carefully. Among estimators built from $N$ independent uniform samples, $1/\sqrt N$ is the rate, and no cleverness changes it. What can change is the constant — that is variance reduction — or the independence, which is what giving up randomness does, and which buys a better exponent at the cost of a much more delicate error analysis.

What is actually being used #

Three settings where the trade decided the method, and none of them is about integrals in the schoolbook sense.

Statistical physics. The average energy of a system of $n$ particles is an integral over a space of dimension $3n$, and $n$ is in the thousands. No grid has ever been contemplated.

Finance. The value of an option depending on a basket of assets over many time steps is an integral whose dimension is assets times steps, routinely in the hundreds.

Rendering. The light arriving at a pixel is an integral over all the paths light could have taken to it, a space with no finite dimension at all. Every ray tracer is a Monte Carlo integrator, and the noise in an under-sampled image is exactly the $\sigma/\sqrt N$ above, made visible.

In all three the alternative is not a slower method; it is no method. That is the sense in which being equally bad in every dimension is a virtue.

What the picture cannot show #

The grid lines are drawn from the exact product formula, so they are smooth; the random line is drawn from measured spreads over twenty-four runs, so it is not. That difference in texture is real information — one is a deterministic quantity and the other is an estimate of a standard deviation — and it reads on the page as one line being better drawn than the other.

The figure also cannot show what happens for an integrand that is not smooth. The grid’s $N^{-2/d}$ assumes two derivatives; for a discontinuous integrand — the indicator of a region, which is what a volume calculation is — the rate collapses to $N^{-1/d}$, and the crossover moves down to two dimensions. So the figure is generous to the grid, and the honest summary is that four is the crossover in the best case for grids and lower otherwise.

And there is nothing here about which random points. The whole calculation assumes independent uniform draws, and a computer supplies neither: it supplies a deterministic sequence chosen to look independent, whose failures have caused real errors in real simulations. Nothing in this figure would detect a bad generator.

An error bar that is itself an estimate #

One further point separates Monte Carlo from a grid, and it is usually counted as an advantage.

A grid gives an answer and no indication of how wrong it is. Bounding the error of a quadrature rule requires bounding a derivative of the integrand, which is normally unavailable and, in the high-dimensional problems where any of this matters, meaningless.

A Monte Carlo run gives an answer and the sample standard deviation of the values it averaged, which estimates $\sigma$, which estimates the error. So the method reports its own uncertainty for free, in the same run, using the same numbers.

The catch is that the estimate is itself random and can be badly wrong in exactly the case where it matters. If the integrand has a rare enormous value — a spike occupying a millionth of the volume and contributing a tenth of the integral — a run of ten thousand points will probably miss it entirely, return an answer that is ten per cent low, and report a small standard deviation, because the sample it saw was genuinely well behaved. The error estimate is only as good as the sample’s coverage of the places that matter, and no amount of arithmetic on the sample can reveal what the sample did not contain.

That failure mode is not hypothetical; it is the standard way Monte Carlo estimates go wrong, and it is the reason importance sampling exists.

The ladder from here #

Rungs above: variance reduction, which attacks the constant rather than the exponent. Importance sampling, stratification and antithetic variates. Markov chain methods, for when the distribution can be evaluated and not sampled. Quasi-random sequences, which beat the square root by abandoning independence. The effective dimension of a problem, which explains why some hundred-dimensional integrals behave like five-dimensional ones. And the error estimate itself, which is a Monte Carlo estimate too and has its own uncertainty.

The shape of the argument #

What decided this comparison was not which method is better but how each one’s cost scales, and the two scalings crossed.

That is the standard form of an asymptotic argument and it is worth noticing how little it says about any particular case. Nothing above establishes that Monte Carlo beats a grid on any specific problem at any specific $N$; it establishes that the exponents differ and therefore that some crossover exists. Where it is depends on constants nobody computed.

The same reasoning appears whenever two costs are compared as a parameter grows, and it comes with the same caveat: a small case can lie about which side of a crossover it is on, and an argument about exponents is silent about everything before the crossing.