Q:

How can I intuit the role of the central limit theorem in breaking the curse of dimensionality for Monte Carlo integration

I would like to more intuitively understand where the power of Monte Carlo integration comes from for large-dimensional domains of integration. 
Other questions on this site have referenced the proof that the scaling of the error in a MC integration goes as $N^{-1/2}$, where $N$ is the number of function evaluations, and does not depend on the dimension $d$ of the domain of integration. On the other hand, the scaling of the error in a uniform sampling integration goes as some power of $N^{-1/d}$. Consequently, for a sufficiently large $d$, one may achieve a desired level of accuracy using a fewer number of function evaluations with a Monte Carlo method than with a uniform numerical quadrature method. The essence of the proof is an invocation of the central limit theorem. I understand the proof on a formal level.
However, I have no intuition for why this proof works. It still seems like "cheating" to me that by either randomly or quasi-randomly selecting the locations to evaluate the function being integrated, you can achieve a more accurate integration than you would by choosing to evaluate the points in a uniform set of spacings (provided that the dimensionality of the integral is large enough). 
I have tried to construct an extremely non-rigorous analogy to conducting a statistical survey. One could either poll people at regular distance spacings in a neighborhood, or instead poll the same number of people in completely randomized locations. If there were no correlation between a person's response and the location where that person lives, then I would be able to conclude that my sampling error would be the same using either method. If there were such a correlation after all, then I might bias my results. Is this sort of reasoning at all on the right track to building intuition for the proof?
To phrase my question another way, it seems that when performing a high-dimensional integral you gain accuracy by including a variety of length scales in the points where you evaluate the functions. Why is that? And, is the role of the randomness essentially to ensure that you use a wide range of length scales?

A:

As you correctly point out, the error in a uniform sampling integration goes as some power of $N^{−1/d}$. The problem of uniform sampling is indeed similar to a non-representative statistical survey. Many of the points in a uniform tensor grid don't add much new information. For example, assume that the function you want to integrate looks like $F(x,y)\approx f(x)+g(y)$ in a local square. If you sample this square with a $k\times k$ tensor grid, you will have used $k^2$ evaluations of $F$, but only received information about $f$ and $g$ for $k$ different arguments. A Taylor expansion of $F$ in a local square as $F(x,y)=ax+by+c+O(x^2+y^2)$ indicates that it is indeed not unusual for $F$ to locally look like that. (Note that no term related to $xy$ occurs in this expression.)
To analyze the connection between uniform tensor grids and the Taylor expansion further, let's look at multivariate polynomials of a fixed degree $k$ in $d$ variables. The function values on a $d$-dimensional uniform tensor grid with $k+1$ points in each direction uniquely determine a polynomial where the max degree of each individual variable (in each term) is smaller or equal than $k$. However, the degree (of an individual term) of a polynomial is defined as the sum of the degrees of the individual variables, so that the determined polynomial has degree $dk$ and $(k+1)^d$ individual terms (=degrees of freedom). On the other hand, a general multivariate polynomial of degree $k$ in $d$ variables only has $\frac{(d+k)!}{d!k!}$ individual terms (=degrees of freedom). For fixed $k$, we have $\frac{(d+k)!}{d!k!}=\frac{d^k}{k!} + O(d^{k-1})$, which is polynomial in $d$ (i.e. not exponential).
What has all this to do with randomness and Monte Carlo integration? It might help explain why the naive tensor product approach doesn't break down for Monte Carlo integration, while it does break down for uniform sampling. Unfortunately, there is no straightforward way to turn an effective (naive) randomized algorithm into an effective (naive) deterministic algorithm. However, we know from experience that problems which can be solved efficiently by randomized algorithms normally can also be solved efficiently by (more complicated) deterministic algorithms. The deterministic high-dimensional integration methods I know work by exploiting recursion (divide and conquer) in a clever way. However, this is not the place to explain hierarchical bases, hierarchical subspace splitting or the Smolyak method. A more black-box approach would be to use space-filling curves as a building block instead, where the recursion is already built into the construction of the curve. I don't know how well this would work for integration, but it has been asked here before. At least I know that the research groups investigating applications of sparse grids often also investigate applications of space filling curves.

