I need to numerically evaluate 2-D integrals of the form: $$\mathcal{I}(\theta) = \int_{0}^{1} \int_0^1 \varphi_\theta(x,y) dx dy$$

where $\varphi_\theta$ is a family of smooth functions indexed by $\theta$; and I need to evaluate $\mathcal{I}(\theta)$ for a large number ($> 10^6$) of $\theta_i \in \Theta$, which are given.

For how the problem is set up, I can evaluate $\varphi_\theta$ on a fine, equispaced 2-D grid (each evaluation involves a few multiplications), and I cannot evaluate $\varphi_\theta$ on other points.

Since the grid is fine, I could simply compute $\mathcal{I}(\theta_i)$ via trapezoidal or Simpson's rule. However, this becomes extremely computationally expensive and unnecessary. For many (but not all) values of $\theta$, $\varphi_\theta(x,y)$ is $\approx 0$ for large part of the integration domain, and/or slowly-varying.

The obvious thing here is to use some sort of adaptive integration method compatible with a fixed equispaced grid. The basic idea I have is to start with a very coarse sub-grid, evaluate the integral and an estimate of the error in each sub-cell, and based on that decide which cells to "zoom in" and recompute with a refined grid. Iterate this a couple of times. There are lots of ways to do this naively, but I am interested in state-of-the-art solutions.

My question are:

• What would be the best (fastest and precise) quadrature approach to compute the (sub)integrals?
• What would be the best method to pick the cells to refine?
• If I use a measure of relative or absolute error, what's the best method to compute that?
• Any alternative idea?

For the record, I am working with MATLAB but I plan to code this part in C via MEX files since I doubt that the adaptive bit can be efficiently vectorized, and I want it to be as fast as possible.

• The "equispaced" requirement makes it difficult, there is no way to drop this? – Bort Aug 26 '16 at 9:34
• @Bort: Not really. – lacerbi Aug 26 '16 at 12:02

I don't think it's quite clear to me why an adaptive scheme would be better. Integration with the trapezoidal rule is very cheap: It's in essence just one addition and one multiplication per grid point. Since you already evaluate your integrand $\varphi_\theta$ at each grid point, the summation to the integral would seem to me to be a completely negligible additional step: surely, evaluation of the integrand at these grid points must be far more expensive than the actual summation.
• Thanks; I clarified this point in the main post. The longer story is that the integrand $\varphi_\theta(x,y)$ is the tensor product and Hadamard product of a number of vectors and matrices which are given on a fine grid. When I say that $\varphi_\theta(x,y)$ is given I mean that I am given the underlying vectors/matrices that then I need to multiply elementwise. So I have a cost for actually evaluating $\varphi_\theta$ (3-4 multiplications). Also, the cost for performing millions of sums ($1000 \times 1000$ grid) is non-negligible when e.g. a $100 \times 100$ grid would work. – lacerbi Aug 25 '16 at 21:06
• But in the case of a tensor product, don't you have $\varphi_\theta(x,y)=a_\theta(x)b_\theta(y)$, which is of course easiest to integrate by two one-dimensional integrals. For the Hadamard (elemntwise) product, the situation is clearly more complicated. – Wolfgang Bangerth Aug 26 '16 at 5:11
• Yep. I have both as a part of the same expression: $\varphi_\theta(x, y) = a_\theta (x) b_\theta(y) c_\theta(x,y)$. Also, I didn't mention that $\theta \in \mathbb{R}^n$, and different (overlapping) subspaces of $\theta$ affect $a$, $b$ and $c$. I would say that these details are irrelevant; the problem I need to solve, in its simplest form, is as stated in my original post. – lacerbi Aug 26 '16 at 5:20