# When is it advantageous to iterate integrals numerically?

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.

• Can you just use different quadrature rules in each dimension? Since $f$ is smooth in $y$, you could use low-order quadrature there and higher-order quadrature along the less-smooth dimensions. The quadrature points and weights be constructed by tensor products of the 1-D rules. That should pretty dramatically cut down the total number of function evaluations. – Tyler Olsen Mar 14 '16 at 14:12
• @TylerOlsen No, but that's the question... (Also, surely high-order rules for smooth dimensions, not the other way?) If the function is sufficiently difficult that I handle the many non-smooth dimensions with QMC, or sparse grids, or tensor products of rules with adaptive refinement, etc., isn't the high precision of integrating over $y$ wasted by the errors in the other dimensions? So why not just use an $(n+1)$-dimensional rule without wasting the evaluations? But I'm not sure what the right general principle is here, so I'm asking if anybody knows something. – Kirill Mar 14 '16 at 16:16
• @TylerOlsen In other words, it should somehow depend on just how easy it is to integrate the function over $y$ (very easy—use $1$-d then $n$-d, but very hard—use $(n+1)$-d). But what's a good guide, or a threshold? – Kirill Mar 14 '16 at 16:20
• Ah, sorry about that. I flipped it around in my head. I can't say that I know a good rule of thumb here, so someone else will have to take this one. – Tyler Olsen Mar 14 '16 at 17:16

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.

• Thank you, but I don't think this answers the question. I've added what I think is a counterexample to the question. Indeed, I think you are almost begging the question when you say things like All of this results in enhanced efficiency for the outer integration routine — it is absolutely not clear that this gains anything for the reason that calculating the integral over $y$ is not free and takes some time itself. So there must be some threshold beyond which it is simply not worth, and the question is where that threshold is. – Kirill Mar 15 '16 at 19:56
• My argument is simply that if you assume an adaptive method with deterministic error control, then evaluating a few dimensions closed-form (or semi closed-form) eliminates inevitable steps that would otherwise be done numerically. But your (excellent) example is one where a deterministic adaptive method would never be used in the first place. – Richard Zhang Mar 15 '16 at 20:19
• Suppose that you did use a standard deterministic adaptive method, e.g. ab-initio.mit.edu/wiki/index.php/Cubature, then I would be very, very surprised if you do not get an entire factor speed-up by cutting out one of the dimensions semi-analytically. – Richard Zhang Mar 15 '16 at 20:25
• Using Cubature exactly how you suggest (which is a perfectly reasonable suggestion) is how I got the "30 times slower" figure in my question in the first place, so I was surprised (hence the question). I only meant the Monte Carlo example as something easy to analyze, I'm actually more interested in deterministic methods. – Kirill Mar 15 '16 at 20:31
• When you say "inevitable steps", I don't think that's really the right way to see this problem. When you integrate out one variable, you oversample greatly in that one variable, out of all the $n+1$ variables, which is something an $(n+1)$-dimensional product rule would never do as it is very inefficient. So it is absolutely not clear that integrating out a variable is always a sensible thing to do. – Kirill Mar 15 '16 at 20:35