R help - Nov 2011 - performance of adaptIntegrate vs. integrate

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?

I'm open to suggestions of alternative ways of dealing with
vector-valued integrands.

Thanks.

baptiste

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
> 
>

The integrand is highly peaked.  It is approximately an impulse function where
much of the mass is concentrated at a very small interval.  Plot the function
and see for yourself.  This is the likely cause of the problem.

Other types of integrands where you could experience problems are: integrands
with singularity at either limit and slowly decaying oscillatory integrands.  As
to why integrate performs better than adaptIntegrate in this situation, I
don't know.  You have to study the two implementations.   Wynn's epsilon
algorithm is an extrapolation method for improving the convergence of a
sequence.  This could be an explanation for the better performance, but I cannot
say for sure.

Hope this is helpful,
Ravi
-------------------------------------------------------
Ravi Varadhan, Ph.D.
Assistant Professor,
Division of Geriatric Medicine and Gerontology School of Medicine Johns Hopkins
University

Ph. (410) 502-2619
email: rvaradhan@jhmi.edu<mailto:rvaradhan@jhmi.edu>


	[[alternative HTML version deleted]]