1
$\begingroup$

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.

Improve this question
$\endgroup$
2
  • $\begingroup$ The "equispaced" requirement makes it difficult, there is no way to drop this? $\endgroup$
    – Bort
    Commented Aug 26, 2016 at 9:34
  • $\begingroup$ @Bort: Not really. $\endgroup$
    – lacerbi
    Commented Aug 26, 2016 at 12:02

1 Answer 1

Reset to default
1
$\begingroup$

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.

It may of course be that you have so many grid points that you need to store the values of the integrands at the grid points on disk, or in a large block of memory, and that that requires loading data back into memory or from memory into the processor cache. The solution in that case is to do the summation in the same loop in which you evaluate the integrand in the first place.

Improve this answer
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$
    – lacerbi
    Commented Aug 25, 2016 at 21:06
  • $\begingroup$ 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. $\endgroup$
    – Wolfgang Bangerth
    Commented Aug 26, 2016 at 5:11
  • $\begingroup$ 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. $\endgroup$
    – lacerbi
    Commented Aug 26, 2016 at 5:20

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.