adaptive Gauss-Kronrod quadrature with vector-valued integrand

Question

So I'm trying to implement a Gauss-Kronrod adaptive quadrature. That is, I want to calculate

$$\int_a^b f(x) dx = \sum_i f(x_i) w_i$$

where f(x) is evaluated at multiple points at once for efficiency. A key point to bear in mind is that $f$ is also vector-valued (possible risk of confusion here with the double meaning, and use of, vectorised); in practice, I'm integrating at once the real part and imaginary part of a complex integrand, or the integrand could be a 3D vector (electric field), etc.

I've worked out the basic quadrature step in R (ultimately to be ported to C++ using the Armadillo library, any Matlab-like matrix syntax will do for the sake of discussion), where I first divide the integral into N portions, and evaluate the weighted sums corresponding to the Gauss and Kronrod nodes.

## nodes and weights borrowed from Octave source code

abscissa = c(-0.9914553711208126e+00, -0.9491079123427585e+00, 
             -0.8648644233597691e+00, -0.7415311855993944e+00, 
             -0.5860872354676911e+00, -0.4058451513773972e+00, 
             -0.2077849550078985e+00,  0.0000000000000000e+00, 
             0.2077849550078985e+00,  0.4058451513773972e+00, 
             0.5860872354676911e+00,  0.7415311855993944e+00, 
             0.8648644233597691e+00,  0.9491079123427585e+00, 
             0.9914553711208126e+00)

weights15 = c(0.2293532201052922e-01,  0.6309209262997855e-01, 
              0.1047900103222502e+00,  0.1406532597155259e+00, 
              0.1690047266392679e+00,  0.1903505780647854e+00, 
              0.2044329400752989e+00,  0.2094821410847278e+00, 
              0.2044329400752989e+00,  0.1903505780647854e+00, 
              0.1690047266392679e+00,  0.1406532597155259e+00, 
              0.1047900103222502e+00,  0.6309209262997855e-01, 
              0.2293532201052922e-01)

weights7  = c(0.1294849661688697e+00,  0.2797053914892767e+00, 
              0.3818300505051889e+00,  0.4179591836734694e+00, 
              0.3818300505051889e+00,  0.2797053914892767e+00, 
              0.1294849661688697e+00)

gauss_kronrod <- function(f, a, b, startN=10, dimf = 2, eps = 1e-3, ...){

  Nquad <- length(abscissa)
  intervals <- seq(a, b, length=startN)
  ## change of variable from [a, b] to [-1, 1]
  shifts <- (intervals[-length(intervals)] + intervals[-1]) / 2
  scalings <- diff(intervals) / 2

  ## scale the weights accordingly
  weightsG <- tcrossprod(weights15, scalings)
  weightsK <- tcrossprod(weights7, scalings)

  ## transformed node positions
  x.scaled <- tcrossprod(abscissa, scalings) +
    matrix(shifts, ncol=length(shifts), nrow=Nquad, byrow=TRUE)

  ## evaluate the function at all points
  fvals <- f(c(x.scaled))

  ## function evaluations are reshaped into a 3D array
  ## rows correspond to dimf values of integrand
  ## columns are the Nquad evalution points in the sub-intervals
  ## slices correspond to the (startN - 1) sub-intervals
  dim(fvals) = c(dimf, Nquad, startN-1)

  ## select which subset of fvals are used for the 7-points kronrod sum
  ikron <- rep(c(FALSE, TRUE), length.out=Nquad)


  gauss <- kronrod <- matrix(0, nrow=dimf, ncol=startN-1)
  for (ii in seq.int(startN-1)){ # integrals for each sub-interval
    gauss[,ii] <- fvals[,,ii] %*% weightsG[,ii]
    kronrod[,ii] <- fvals[,ikron,ii] %*% weightsK[,ii]
  }
  ## relative difference between the two rules 
  errors <- abs(gauss - kronrod) / abs(gauss)
  test <- which(apply(errors, 2, max) > eps)

  ## net integral, just sum the partial sums...
  rowSums(gauss)

}

# a vector-valued integrand...
f <-  function(x) rbind(exp(x/(x+2))*sin(x), cos(x))

gauss_kronrod(f, 2, 110, startN=10)

My problem comes when trying to implement an adaptive scheme; basically, at each iteration I would like to compare the values of gauss vs kronrod for each of the subintervals, and from there decide where to subdivide further. What would be a good strategy to store these intermediate results during the iteration, and keep track of which subintervals sums go where etc. in the final summation?

More general advice on this endeavour will also be appreciated, I'm not usually keen to implement such generic operations myself, but it seems it might be easier than linking different libraries together (I was contemplating cubature, initially, because my end goal is 2D integration).

Answer 1

I suggest to read the Quadpack book (Quadpack, Piessens, Robert; de Doncker-Kapenga, Elise; Überhuber, Christoph W.; Kahaner, David (1983). QUADPACK: A subroutine package for automatic integration. Springer-Verlag. ISBN 978-3-540-12553-2) and Pedro Gonnet's PhD thesis "Adaptive Quadrature Re-revisited (available as pdf here). Pedro is a contributor to SciComp SE.

What typically is done is a priority queue data structure:

1. You apply Gauss-Kronrod on the whole initial interval and estimate the error
2. You split the interval symmetrically in two and again apply Gauss-Kronrod on the intervals, estimate the errors and add a pair (or tuple) "interval, error" to the priority queue
3. If the total error is smaller than the requested tolerance (or if you reached the maximal amount of subdivisions), you stop
4. If not, you pop the top element from the queue (highest error), you split it in two and apply Gauss-Kronrod, add those two new elements to the priority queue and go back to 3.

Answer 2

Well, I've thought of a strategy, which is the following. First, I create a function to take care of calculating the partial sums over a list of intervals, and checking which are good enough to keep. Second, for those intervals that were not good, a helper function splits them in two. Finally, at each step the good contributions are added to the final result. Here's an implementation in R,

gk_segments <- function(f, segments=list(c(2, 50), c(50, 100)),  dimf = 2, 
                        eps = 1e-3, ...){

  Nquad <- length(abscissa)
  Nsegments <- length(segments)
  shifts <- sapply(segments, sum) / 2
  scalings <- sapply(segments, diff) / 2
  weightsG <- tcrossprod(weights15, scalings)
  weightsK <- tcrossprod(weights7, scalings)

  x.scaled <- tcrossprod(abscissa, scalings) +
    matrix(shifts, ncol=length(shifts), nrow=Nquad, byrow=TRUE)

  fvals <- f(c(x.scaled), ...)

  ## function evaluations are reshaped into a 3D array
  ## rows correspond to dimf values of integrand
  ## columns are the Nquad evalution points in the sub-intervals
  ## slices correspond to the Nsegments sub-intervals
  dim(fvals) = c(dimf, Nquad, Nsegments)

  ikron <- rep(c(FALSE, TRUE), length.out=Nquad)

  gauss <- kronrod <- matrix(0, nrow=dimf, ncol=Nsegments)
  for (ii in seq.int(Nsegments)){ # integrals for each sub-interval
    gauss[,ii] <- fvals[,,ii] %*% weightsG[,ii]
    kronrod[,ii] <- fvals[,ikron,ii] %*% weightsK[,ii]
  }

  errors <- abs(gauss - kronrod) / abs(gauss)
  test <- which(apply(errors, 2, max) > eps)

  if(!length(test))
    return( list(converged = rowSums(gauss),
                 todo = list()))

  list(converged = rowSums(gauss[ ,-test]),
       todo = segments[test])
}

# f <-  function(x) rbind(exp(x/(x+2))*sin(x), cos(x))
# gk_segments(f, segments=split_interval(c(2, 110), 10))

# split a given interval [x0, x1] into a list of N subintervals
split_interval <- function(x, N=2){

  breaks <- seq(x[1], x[2], length=N+1)
  lapply(seq.int(N), function(ii) c(breaks[ii], breaks[ii+1]))
}

gk_adaptive <- function(fun, a, b, startN=10, dimf = 2, 
                         eps = 1e-3, maxiter = 100, ...){

  ## first step
  segments <- split_interval(c(a, b), startN)
  tmp <- gk_segments(f=fun, segments=segments, dimf=dimf, eps=eps, ...)

  result <- tmp$converged # partial sums that have already converged
      segments <- unlist(lapply(tmp$todo, split_interval, N=2), recursive=FALSE)
  iter <- 1

  ## we stop if there's no subinterval left to split
  ## or we've reached the maximum number of iterations
  finished <- (!length(tmp$todo)) || (iter > maxiter)

  while(!finished){
    message("iteration #", iter)
    ## work on the new subsegments
    tmp <- gk_segments(f=fun, segments=segments, dimf=dimf, eps=eps, ...)
    result <- result + tmp$converged ## add the good chunks
        segments <- unlist(lapply(tmp$todo, split_interval, N=2), recursive=FALSE)

    iter <- iter + 1
    stopping <- iter > maxiter
    if(stopping)
      message("Maximum number of iteration reached, result is likely wrong")
    finished <- (!length(tmp$todo)) || stopping
  }

    return(result)
}

f <-  function(x) rbind(sin(x), cos(x))
gk_adaptive(f, 2, 110)
gk_adaptive(f, 2, 110, startN=10, eps=1e-10)
gk_adaptive(f, 2, 110, startN=3, eps=1e-10, maxiter=2)

I don't think it will be easy to generalise to 2 or more dimensions, because choosing the right square / cubes to focus on iteratively seems harder to code. As a less efficient alternative, I can just chain multiple integrations in 1D, keeping all but one of the variables constant at a time.