Convergence of Monte Carlo integration

Question

In my research, one of the steps is to choose a numerical method to estimate $\int_a^b f(t)dt$, where $f$ is Lipschitz continuous but not differentiable. For simplicity, I used midpoint rule but the overall result is not that good. I was thinking whether the problem came from the fact that I didn't select a correct numerical quadrature rule.

Based on my understanding, to guarantee convergence, most deterministic quadrature rules require the integrand be certain order differentiable. Since in this case, $f$ is not smooth enough, it's reasonable that the mid-point rule can't guarantee convergence.

My supervisor suggested me to use Monte-Carlo integration. But I am not sure about the convergence condition for this numerical method. I only know regardless of the dimension, this method has a $\frac{1}{\sqrt{N}}$ convergence order. Does this result hold for all $L^1$ function? Do we need any additional assumption for integrand $f$ for the convergence? Could you provide me some resources about the convergence analysis for Monte Carlo integration?

Answer 1

Let the MC integral estimate be $$ S_n = \frac{b-a}{n}\sum_{1\leq k\leq n} f(x_k), $$ where $x_k$ are i.u.d. on the interval $[a,b]$. So long as the function $f$ is $L^1$, the mean exists, and $$\mathbb{E}[S_n] = I, \qquad I = \int_a^b f(x)\,dx, $$ so by the law of large numbers $S_n$ will converge to $I$.

To get the $1/\sqrt{n}$ convergence, the function also needs to be square-integrable: $$ \mathbb{V}[S_n] = \frac{|b-a|}{n}\int_a^b \big(f(x)-\bar f\big)^2\,dx, $$ (where $\bar f$ is the function average on $[a,b]$), so if $f$ is also $L^2$, then the central limit theorem applies, and $S_n$ would be normally distributed with mean $I$ and variance $n^{-1}|b-a|\mathbb{V}[f]$. When the function is not square integrable, the law of large numbers still applies, but this error distribution is no longer valid, it would have a fat tail.

Lastly, non-differentiable functions are not that bad for numerical integration, so long as the singularities are easy to locate, and the algorithm performs some kind of an adaptive interval bisection procedure that isolates the singularities. Alternatively, a common trick is to split the interval by hand so that all singularities are at subinterval endpoints, and invoke quadrature schemes on individual subintervals. I believe most mature quadrature libraries can handle reasonable singularities quite well. In particular, one of the best quadrature schemes, the double-exponential rule, has no problems at all with regular singularities at the interval endpoints. Others, like Gauss-Legendre rules, usually have modified forms to address singularities.

You don't say what your function $f$ is, but my opinion would be that it has to be pretty pathological (e.g., not differentiable anywhere) before MC methods would outperform a mature quadrature scheme in 1d. Not midpoint, though — you'd never expect that to work well.