When is it advantageous to iterate integrals numerically?

Question

If there is an $(n+1)$-dimensional integral of the form $$ \int_{[0,1]^{n+1}} f(x, y)\,\mathrm{d}^n x \,\mathrm{d}y,$$ normally one would evaluate this using a multi-dimensional integration library over the whole domain, $[0,1]^{n+1}$.

But are there some conditions where it might make sense to perform the integral over $y$ separately, using a one-dimensional quadrature, and then use the multi-dimensional integration library to evaluate the integrand over the other $n$ coordinates? $$ \int_{[0,1]^n}g(x)\,\mathrm{d}^nx, \qquad g(x) = \int_0^1 f(x,y)\,\mathrm{d}y. $$

This might make sense, for example, if $f$ is especially smooth as a function of $y$, but not $x$. But how smooth exactly does it have to be in this case? My guess was that it almost never makes sense because far too many of the 1-d quadrature evaluation points would be "wasted", but I'm not so sure this always applies. Is this guaranteed by the design of the high-dimensional integration methods?

In my own case, $f$ is black-box, but piecewise-smooth in $y$, and has an unknown amount of kinks and jumps in $x$ at unknown locations, and $n$ is quite high ($n\geq 4$), so the integral over $x$ has to be done with something specifically for many dimensions. The integral over $y$ can be done with something regular like quadgk. In this example, the function is smooth enough in $y$ that it almost seems to work, but repeated integration ends up being 30 times slower, so I am wondering if the approach is misguided.

If you know where this is already discussed in the literature, that would be helpful too.

Example. (of why this is not trivial) Consider an "easy" integral, which is very smooth, unlike what I'm really interested in: $$ \int_{[0,1]^n} e^{-x_1 x_2\cdots x_n}\,\mathrm{d}^n x = \mathrm{F}\left(\begin{array}{}\{1,\ldots,1\}_n\\\{2,\ldots,2\}_n\end{array}\middle| -1\right). $$ We could do naive $n$-dimensional Monte Carlo on the integrand, or naive $(n-1)$-dimensional Monte Carlo on the integrand integrated once w.r.t. $x_1$, which is $g(x_{2:n}) = -(e^{-a}-1)/a$ (where $a=x_2\cdots x_n$).

With some algebra, I calculated that the variance of $(n=5)$-dimensional $N$-point MC estimate is $0.00244N^{-1}$, and it is $0.00167N^{-1}$ for the $4$-dimensional integral of $g$, for a variance reduction by a factor of $1.5$.

This is a paltry reduction in variance: it would be negated by using $1.5$ times as many samples points, and this is offset by the fact that the internal integrand might be more than $1.5$ times slower to evaluate. If the function $g = (1-e^{-a})/a$ above happens to be more than $1.5$ times slower, this represents a net loss of accuracy, keeping computation time fixed.

Presumably the same kind of tradeoff applies when considering a deterministic rule for integrating over $x$. The Monte Carlo method makes this analysis much easier than the general case, because integrating over $y$ acts like a very straightforward variance reduction technique. But I'm really much more interested in deterministic methods, which I haven't been able to analyze as easily.

Answer 1

Clarification: My answer is specifically written for adaptive integration routines with deterministic error control like this one. It becomes moot for sparse grid and Monte Carlo based integration routines, whose error control is not performed in the manner described below.

A significant cost of automatic black-box, tensor-product-based integration routines is error control, from two aspects

1. Wasted function evaluations. All adaptive integration work by estimating the integrand and the error using a low-order rule or coarser partitioning, and repeating this with progressively higher-order rules or finer partitions until the error requirements are satisfied. Nested integration rules allow some of work done in previous steps to be recycled, but often not all.
2. In an effort to conserve function evaluations, highly nested rules like Gauss-Kronrod or Newton-Cotes are often used in adaptive integration codes. Nested quadrature rules are suboptimal quadrature rules, in that they are considerably less accurate than optimal rules (e.g. Gauss-Legendre and Clenshaw Curtis) for a particular class of functions over a fixed quadrature order. In other words, nested rules make less efficient use of function evaluations.

The inefficiencies associated with error control in the $y$ dimension can be entirely removed using a semi-analytic treatment over $y$. Suppose that $f(x,y)$ were smooth (i.e. not just piecewise smooth) in $y$. We are free to pick an optimal quadrature rule, and to derive an a prior bound on the quadrature order $r$ in order to guarantee $$\left|\sum_{k=1}^r w_k f(x,y_k) - \int_{[0,1]}f(x,y)\,\mathrm{d}y\right| \le \epsilon\qquad \text{for all }x\in[0,1]^n,$$ using, say, the Lipschitz constant of $f(x,y)$ over $y$ for all $x$. In a semi-analytical treatment of variable error control, we may even let the quadrature order $r$ vary as a function with respect to $x$ and $\epsilon$ to further enhance efficiency.

Replacing numerical steps by optimal, semi-analytical ones results in enhanced efficiency for the outer integration routine. Considering that $g(x)$ lies at the heart of an $n$-dimensional integral, any small gains made in the evaluation of $g(x)$ will be greatly magnified over the full integration.

In the case where $f(x,y)$ is only piecewise smooth, optimal quadrature rules can be developed so long as the boundaries of the piecewise segments are known, so similar arguments hold. Even if the boundaries are not known, heuristically, Clenshaw-Curtis is often able to handle piecewise smooth functions very well.

To give an application example, this exact issue arose for me in the evaluation of volume-to-volume singular integrals in this paper, and my treatment is similar to the one proposed above. As a rule of thumb, it is always advisable to remove as many dimensions as possible using analytical arguments before feeding the problem through a black-box integration routine.