Implementing AbsTol for quadcc

Previous Topic Next Topic
 
classic Classic list List threaded Threaded
4 messages Options
Reply | Threaded
Open this post in threaded view
|

Implementing AbsTol for quadcc

Rik-4
9/22/17

I've been working with Nicholas on bug #52074
(https://savannah.gnu.org/bugs/?52074) to implement Matlab-compatible
functions integral2, integral3 which in Octave are wrapper functions which
call down to one of the quadrature routines quad, quadl, quadv, quadgk, or
quadcc.

The default integrator in Octave is quadcc which is a very good all-purpose
quadrature routine in that it can handle discontinuities, singularities,
NaNs, Infinite lower or upper bounds, etc.  It also happens to be written
in C++ so it is faster than quadl, quadv, and quadgk.  However, the
algorithm for quadcc only implements relative tolerance RelTol.  The
calling form is quadcc (f, a, b, reltol).  Because quadcc is not a Matlab
function, there is no requirement to have the same functional interface,
but all of the other routines use quadXX (f, a, b, abstol).  If you want to
swap back and forth between different integrators because of certain
properties of your integrand then this is inconvenient.  It also breaks
good UX conventions which favor consistency.

In addition, by only implementing RelTol there are some integrals which
become pathological.  For example, bug report #46439
(https://savannah.gnu.org/bugs/?func=detailitem&item_id=46349) shows that a
periodic function, integrated over exactly one complete period, can take an
exhaustively long time.  This code:

tic; [q] = dblquad(@(x,y) sin(2 * pi * x ) * sin(2 * pi * y),0,1,0,1, [],
'quadcc'); toc

takes 35 seconds, versus .163 seconds for quadgk.

It wasn't hard, so I made a patch to implement both absolute and relative
tolerance.  Now, the double integral above completes in .013 seconds or 12X
faster than quadgk.  I wanted to have access to both tolerances for the
algorithm so I did this

+The optional argument @var{tol} is a 1- or 2-element vector that specifies the
+desired accuracy of the result.  The first element of the vector is the
desired
+absolute tolerance, and the second element is the desired relative tolerance.
+To choose a relative test only, set the absolute tolerance to zero.  To choose
+an absolute test only, set the relative tolerance to zero.  The default
+absolute tolerance is 0 and the default relative tolerance is @math{1e^{-6}}.

This is similar to the interface for quad().

How do people feel about making this change?  Should I default it to using
[RelTol AbsTol] so existing code can run unmodified?  I'd prefer not to
since the interface would then differ from the other quad routines.  Should
there be a warning in place for two versions of Octave when the single
tolerance syntax is used informing the user that this value will now be
used for AbsTol rather than RelTol?  Should the AbsTol default to something
more sensible than 0, such as 1e-10 which is what quadgk uses?

--Rik

Reply | Threaded
Open this post in threaded view
|

Re: Implementing AbsTol for quadcc

NJank
On Fri, Sep 22, 2017 at 1:31 PM, Rik <[hidden email]> wrote:
9/22/17

I've been working with Nicholas on bug #52074
(https://savannah.gnu.org/bugs/?52074) to implement Matlab-compatible
functions integral2, integral3 which in Octave are wrapper functions which
call down to one of the quadrature routines quad, quadl, quadv, quadgk, or
quadcc.


To elaborate: 

Matlab has apparently deprecated most of it's quad style integrators. they released integral (and 2 and 3) in 2012, and looking through the pages for the quad-style functions, most say "will be removed in a future release. Use integral instead.".Applies to:

quad, quadl, quadv, dblquad, and triplequad (but not quadgk or quad2d)

Note the two not deprecated and integral and 2 and 3 all make use of the same general calling form, like we use now for quadgk:

quadgk(f, lim1, lim2, param1, val1, param2, val2, ...)

In addition to some other parameters, all are expected to accept both abstol and reltol, which are used in the same way specified for quadgk. Trying to have the wrappers handle abstol & reltol compatibly, with appropriate defaults, and having the underlying integrators actually respect those inputs has made the wrappers a bit convoluted.

The need to call both abstol and reltol limits integral to mostly using quadgk (with quadv for array valued functions). dblquad and triplequad only pass a single tol (incorrectly described in doc as abstol since it defaults to quadcc) to the quad specified, and they don't support passing any additional params to quad.

having quadcc accept a single vector parameter would make integral2 and 3 a bit faster and easier, and almost fully compatible. I am concerned about the backwards compatibility.

Full compatibility on integral will likely require someone to recode an arrayvalued version of quadgk or quadcc. For integral2 and 3, we need a quad2d equivalent.  But this would be a step in the right direction.

Nick J.
Reply | Threaded
Open this post in threaded view
|

Re: Implementing AbsTol for quadcc

Marco Caliari-4
In reply to this post by Rik-4
On Fri, 22 Sep 2017 10:31:39 -0700
Rik <[hidden email]> wrote:

> How do people feel about making this change?  Should I default it to using
> [RelTol AbsTol] so existing code can run unmodified?  I'd prefer not to
> since the interface would then differ from the other quad routines.  Should
> there be a warning in place for two versions of Octave when the single
> tolerance syntax is used informing the user that this value will now be
> used for AbsTol rather than RelTol?

Yes, maybe with the hint on how to disable that specific warning.

> Should the AbsTol default to something
> more sensible than 0, such as 1e-10 which is what quadgk uses?

Yes, this is the default in integral?, too.

Marco

Reply | Threaded
Open this post in threaded view
|

Re: Implementing AbsTol for quadcc

NJank
Yes, this is the default in integral?, too.


Yes, matlabs default tolerances all now seem to be abstol = 1e-10, reltol = 1e-6. This matches what octave uses for quadgk, with the exception that for 'single' types octave uses 1e-5 and 1e-4. Matlab keeps the higher tol for singles and the online documentation just warns that you may need to use a larger value for convergence. The thought was that octave would apply its current behavior with all the other functions, making sure it's detailed in all the docs, maybe highlighted with a 'known incompatibility' note.