[R] performance of adaptIntegrate vs. integrate
Hans W Borchers
hwborchers at googlemail.com
Fri Nov 11 10:44:57 CET 2011
baptiste auguie <baptiste.auguie <at> googlemail.com> writes:
>
> Dear list,
>
> [cross-posting from Stack Overflow where this question has remained
> unanswered for two weeks]
>
> I'd like to perform a numerical integration in one dimension,
>
> I = int_a^b f(x) dx
>
> where the integrand f: x in IR -> f(x) in IR^p is vector-valued.
> integrate() only allows scalar integrands, thus I would need to call
> it many (p=200 typically) times, which sounds suboptimal. The cubature
> package seems well suited, as illustrated below,
>
> library(cubature)
> Nmax <- 1e3
> tolerance <- 1e-4
> integrand <- function(x, a=0.01) c(exp(-x^2/a^2), cos(x))
> adaptIntegrate(integrand, -1, 1, tolerance, 2, max=Nmax)$integral
> [1] 0.01772454 1.68294197
>
> However, adaptIntegrate appears to perform quite poorly when compared
> to integrate. Consider the following example (one-dimensional
> integrand),
>
> library(cubature)
> integrand <- function(x, a=0.01) exp(-x^2/a^2)*cos(x)
> Nmax <- 1e3
> tolerance <- 1e-4
>
> # using cubature's adaptIntegrate
> time1 <- system.time(replicate(1e3, {
> a <<- adaptIntegrate(integrand, -1, 1, tolerance, 1, max=Nmax)
> }) )
>
> # using integrate
> time2 <- system.time(replicate(1e3, {
> b <<- integrate(integrand, -1, 1, rel.tol=tolerance, subdivisions=Nmax)
> }) )
>
> time1
> user system elapsed
> 2.398 0.004 2.403
> time2
> user system elapsed
> 0.204 0.004 0.208
>
> a$integral
> > [1] 0.0177241
> b$value
> > [1] 0.0177241
>
> a$functionEvaluations
> > [1] 345
> b$subdivisions
> > [1] 10
>
> Somehow, adaptIntegrate was using many more function evaluations for a
> similar precision. Both methods apparently use Gauss-Kronrod
> quadrature, though ?integrate adds a "Wynn's Epsilon algorithm". Could
> that explain the large timing difference?
Cubature is astonishingly slow here though it was robust and accurate in
most cases I used it. You may write to the maintainer for more information.
Even a pure R implementation of adaptive Gauss-Kronrod as in pracma::quadgk
is faster than cubature:
library(pracma)
time3 <- system.time(replicate(1e3,
{ d <<- quadgk(integrand, -1, 1) } # 0.0177240954011303
) )
# user system elapsed
# 0.899 0.002 0.897
If your functions are somehow similar and you are willing to dispense with
the adaptive part, then compute Gaussian quadrature nodes xi and weights wi
for the fixed interval and compute sum(wi * fj(xi)) for each function fj.
This avoids recomputing nodes and weights for every call or function.
Package 'gaussquad' provides such nodes and weights. Or copy the relevant
calculations from quadgk (the usual (7,15)-rule).
Regards, Hans Werner
> I'm open to suggestions of alternative ways of dealing with
> vector-valued integrands.
>
> Thanks.
>
> baptiste
>
>
More information about the R-help
mailing list