bugGNU Octave - Bugs: bug #62412, quadcc does not warn when...

 
 

bug #62412: quadcc does not warn when quadrature fails to converge

Submitter:  None
Submitted:  Wed 04 May 2022 02:31:12 PM UTC
   
 
Category:  Octave Function Severity:  2 - Minor
Priority:  5 - Normal Item Group:  Missed Error or Warning
Status:  Fixed Assigned to:  None
Originator Name:  Originator Email:  -email is unavailable-
Open/Closed:  * Closed Release:  * dev
Operating System:  * Any Fixed Release:  None
Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Sun 15 May 2022 07:44:03 PM UTC, comment #27: 

The remaining issue of "ArrayValued" integration is now the subject of bug #62468.

Marking this bug as Fixed and closing report.

Rik <rik5>
Group administrator
Wed 11 May 2022 10:31:25 PM UTC, comment #26: 

I think the more useful documentation is for integral at https://www.mathworks.com/help/matlab/ref/integral.html.

This only allows one scalar for absolute tolerance and another for relative tolerance.  There is no change when the property "ArrayValued" is true so presumably the same tolerances are applied to every column of the output.

Rik <rik5>
Group administrator
Wed 11 May 2022 09:53:30 PM UTC, comment #25: 

matlab's public facing documentation for quadv is here:
https://www.mathworks.com/help/matlab/ref/quadv.html

that's really all we have to go off of. no cited references either, unless older versions of the matlab documentation can be found that provided more detail.

Nicholas Jankowski <nrjank>
Group Member
Wed 11 May 2022 09:20:17 PM UTC, comment #24: 

Yes, I think I can do that. Can you please point me to Matlab's specification of vector-valued integration? A question is for instance how the tolerances are specified -- if the elements of the vector-valued function to be integrated have different characteristic scales, it would be natural to specify a vector of absolute tolerances.

Michael Leitner <mleitner>
Wed 11 May 2022 07:02:06 PM UTC, comment #23: 

I checked in a change to have quadgk behave in the same way as quadcc and silence the warning when a second output is requested (https://hg.savannah.gnu.org/hgweb/octave/rev/1fcfe9443917).

What remains seems to be how to have vectorized, "ArrayValued", integrations.  I don't think we should modify quadv.m to have an ERR output.  For Matlab compatibility, the second output must be the number of function evaluations.  In addition, the Matlab documentation says that quadv() will be removed in a future version and recommends using integral() with the "ArrayValued" option.

Michael speculated that it would not be too difficult to improve quadgk.m to support this feature.  I think that is the best course and should become a new issue report.

Rik <rik5>
Group administrator
Tue 10 May 2022 07:08:44 PM UTC, comment #22: 

thanks for the cleanup. my opinions on the outstanding questions:

1 - I think consistent behavior makes sense, especially if they may be called by integral.

2 - i think if the underlying integrators are consistent, the warnings should be left there, with the notion that they're best handled by the routine doing the calculation. with the warning 'living in' the integrator, it will also tell the user exactly which function is being used/issuing the warning vs it just giving a 'warning: integral: ..."

3 - yes. I assume for matlab and backward compatability it will have to be done with an optional 3rd argument. but that will be an easy enough change to integral. As Michael pointed out in comment #12, though, the actual error calc may take a bit of work to do right.

Nicholas Jankowski <nrjank>
Group Member
Tue 10 May 2022 06:45:38 PM UTC, comment #21: 

Forgot to document this, but I reviewed your patch, made a few changes for Octave coding conventions, and checked it in here: https://hg.savannah.gnu.org/hgweb/octave/rev/347fd0653508.

The issue in the title of the bug report is now resolved.

Further questions which were discussed, but I didn't see a definitive answer on are:

1) Should quadgk be changed to not emit a warning about convergence if the second output ERR is requested?  This would make it consistent with quadcc.  Nicholas already has a patch for this.

2) Should integral.m be changed to always call the underlying quadrature routine with 2 outputs and then issue a consistently worded warning message?  I will note that the language from the quadgk warning was copied into the quadcc message so they are consistent right now already.

3) Should a new feature request be filed on Savannah to modify quadv.m to return an estimate of the error so that it would be consistent with other quadrature routines?

Rik <rik5>
Group administrator
Mon 09 May 2022 03:15:29 PM UTC, comment #20: 

ok, thanks for the pointers.  I believe the attached patch (untested) should create a conditional error message that displays if the user does not request err as an output.

(file #53183)

Nicholas Jankowski <nrjank>
Group Member
Sat 07 May 2022 04:50:35 PM UTC, comment #19: 

Alas, it is true that some of this is not documented well.  For an example, see addpath in load-path.cc.  The DEFMETHOD declaration is


DEFMETHOD (addpath, interp, args, nargout,


and later in the code it is used as


  if (nargout > 0)
    retval = lp.path ();


"nargout" is the standard name to use for the 4th argument so that the C++ code resembles what m-file code would look like.

Rik <rik5>
Group administrator
Sat 07 May 2022 11:46:36 AM UTC, comment #18: 

The fourth argument of the DEFMETHOD macro can be used as the number of output arguments in the definition of built-in functions.
See: https://hg.savannah.gnu.org/hgweb/octave/file/760d38460452/libinterp/corefcn/defun.h#l111

#define DEFMETHOD(name, interp_name, args_name, nargout_name, doc)      \
  DECLARE_METHOD (name, interp_name, args_name, nargout_name)

//! Macro to define a builtin method with certain internal name.
//!
//! @warning Consider to use #DEFMETHOD, unless you have good reason.
//!
//! For detailed information, see \ref Macros.
//!
//! This macro can be used when @p name cannot be used directly (for example if
//! it is already defined as a macro).  In that case, @p name is already a
//! quoted string (thus unaffected by macros), and the internal name of the
//! method is given by @p fname.
//!
//! @param name The **quoted** name of the method that should be callable
//!             by the interpreter.
//! @param fname The internal **unquoted** name of the method.  This internal
//!              name is by convention prepended by an 'F'.
//! @param interp_name The name of the 'octave::interpreter' reference that can
//!                    be used by this method.  If this value is omitted,
//!                    there is no access to the interpreter and one should
//!                    use #DEFUNX to define a function instead.
//! @param args_name The name of the octave_value_list variable used to pass
//!                  the argument list to this method.  If this value is
//!                  omitted, the method cannot access the argument list.
//! @param nargout_name The name of the 'int' variable used to pass the number
//!                     of output arguments this method is expected to
//!                     produce from the caller.  If this value is
//!                     omitted, the method cannot access this number.
//! @param doc Texinfo help text (docstring) for the method.
//!
//! @see DEFMETHOD


(That info is probably also available somewhere else in a nicer format. But I don't know where. 🙂)

Markus Mützel <mmuetzel>
Group administrator
Fri 06 May 2022 07:00:37 PM UTC, comment #17: 

quadcc should just require the same nargout check, but not sure how that's done in C.  I see nargin gets defined from args.length. what's the nargout equivalent?

Nicholas Jankowski <nrjank>
Group Member
Thu 05 May 2022 08:42:41 PM UTC, comment #16: 

ok. if that's the case than for quadgk then just changing the final warning to:


  if ((nargout > 1) && (err > max (abstol, reltol * abs (q)))))
    warning (warn_id,
             "quadgk: Error tolerance not met.  Estimated error %g", err);
  endif


simple patch for that attached.

there are other warnings, but they all terminate iteration as well. should to stay even when the end tolerance warning is turned off? Includes:
"quadgk: maximum interval count (%d) exceeded"
and
"quadgk: non-finite integrand encountered")


(file #53177)

Nicholas Jankowski <nrjank>
Group Member
Thu 05 May 2022 07:32:40 PM UTC, comment #15: 

Yes. If an output is requested that would indicate a failure, then there is no need for an additional message. However, in a strict sense only quad (with its IER output) does that. quadgk and quadcc give the expected error, but the user would have to compare that to the requested accuracy. Would that count according to your criterion "an output that would allow to determine..."?

Michael Leitner <mleitner>
Thu 05 May 2022 07:26:29 PM UTC, comment #14: 


>  it doesn't check if nargin>1 first, but i noticed quadgk warns whether err is requested or not, so figured this could do similarly.


Actually, that was behavior I wanted to change in quadgk.

See this note from the NEWS file from version 7.1


- Many functions in Octave can be called in a command form---no
parentheses for invocation and no return argument assignment---or in a
functional form---parentheses and '=' for assignment of return values.

    **Command Form Example**

    `mkdir new_directory`

    **Function Form Example**

    `status = mkdir ("new_directory")`

    Octave now handles errors that occur in a consistent manner.  If
    called in command form and there is a failure, an error is thrown
    and a message printed.  If called in functional form, no error or
    message is printed and the failure is communicated to the programmer
    via the output status variable.


I would like to have the same consistency here.  If the user hasn't asked for an output that would allow them to determine that the algorithm failed then they need to be given a warning.  Otherwise, the warning should be suppressed and the programmer will handle interpreting the results and what should be done (issue warning?  issue error?  something else?).

Rik <rik5>
Group administrator
Thu 05 May 2022 05:21:25 PM UTC, comment #13: 

Re: adding a quadcc warning, I think just the following should suffice just before the final return statements:


  if (err > std::max (abstol, fabs (igral) * reltol)
    warning ("quadcc: Error tolerance not met. Estimated error: %e\n", err);


the attached patch does that, but since i'm not set up to build and test didn't want to push it into quadcc.  it doesn't check if nargin>1 first, but i noticed quadgk warns whether err is requested or not, so figured this could do similarly.

(file #53176)

Nicholas Jankowski <nrjank>
Group Member
Thu 05 May 2022 03:01:58 PM UTC, comment #12: 

I think if it is only a question of complex integrands, then any integrator should do: just call it once with the real part of the integrand and once with the imaginary part, and combine those two afterwards. In fact, one could argue that one should do this always, even when calling the integrators themselves: why should it not be possible to integrate, e.g., sqrt(x^2-1) from -2 to 2? Of course, it is different if you are talking about integration over paths in the complex plane.

I looked into quadv.m, and unless I am strongly mistaken, it computes its internal error incorrectly (or put more positively, with wildly optimistic assumptions), see the FIXME in line 193. So in general it doesn't even follow its tolerance target. And with respect to its ability to handle vector-valued integrands: I think it would not be hard at all to modify quadgk to allow vector-valued integrands.

It seems that in quadgk the errors are really just the sums of the absolute differences in the Gauss and Kronrod steps of the accepted subintervals. Thus, it can often be the case that this is a very conservative estimate if there is cancellation.


Michael Leitner <mleitner>
Thu 05 May 2022 02:06:39 PM UTC, comment #11: 

"It would be nice if in addition, when called from integral, the different called quadrature functions could be muted, and the estimated returned error compared within integral to return a warning that is formatted agnostically to which integrator is used, but this is not necessary."

the 7.1.0 release version now routes complex integrands to quadgk after bug #58636, which supposedly can intelligently(?) handle them. Do you think there may be a better way to process that?

I agree that a single , uniform error estimate would be best.  quadv currently doesn't have a typical estimate as an output, though. just an evaluation count, which can sort-of be used as an arbirtrarily scaled convergence metric. I'm sure it has an internal error metric, but I haven't looked that closely.  It would be just as easy to add an optional third output with that metric without breaking compatability if it was available as a simple scalar. Since it was the only integrator we had that could handle vector-valued functions, its inclusion was needed to maximize the wrapper's compatability.

for better algorithm documentation:

quadcc: states it provides "an estimate of the absolute integration error", and doesn't describe how it handles sigularities or inf limits but does provide the reference for Gonnet where a better definition could be extracted if its not obvious from the code (after verifying we correctly followed the method).
quadgk: describes error as "ERR is an approximate bound on the error in the integral 'abs (Q - I)', where I is the exact value of the integral.", but it could probably better describe how that's calculated. Also no informatino given on transformations but reference to Shampine is given.
quadv: neither is given, nor a reference, but it says it's Simpsons with  recursive subdivision.

Of course, Octave's integral wrapper itself was really meant to just be a bandaid pending an actual integral rewrite. From the matlab help for integral, "quadgk and integral use essentially the same integration method" it looks like they just expanded quadgk to add array valued function handling, as almost all other options are the same. Maybe we over complicated things trying to keep quadcc for performance.

Nicholas Jankowski <nrjank>
Group Member
Thu 05 May 2022 08:52:00 AM UTC, comment #10: 

Some points:

Yes, it is a good idea to make quadcc behave like the others and check for convergence. It would be nice if in addition, when called from integral, the different called quadrature functions could be muted, and the estimated returned error compared within integral to return a warning that is formatted agnostically to which integrator is used, but this is not necessary.

Handling complex integrands would obviously be no problem, within integral you could pass in the real and imaginary parts of the integrand separately (which would be a good idea anyway, as this could lead to better efficiency of the adaptive refinement, when the hard parts of real and imaginary components are different).

In principle, a main point to consider would be the fact that Gaussian quadrature (done by quadgk) and Clenshaw-Curtis quadrature (quadcc) have different behaviour. However, surprisingly in typical cases they do seem to have about equal efficiency: https://doi.org/10.1137/060659831. Thus, I think that also the present overview in the documentation is not really correct: perhaps this is about the default accuracy, but in principle, as long as the integrand is not too obstinate, it should be possible to get close to eps with all these algorithms. The only significant distinction, I would say, is between quadv and the rest: the former uses the Simpson rule, which for smooth enough integrands has much lower order (needs much more function evaluations for high accuracy) than all the other ones.

The opinions in the paragraph above are with respect to the characteristics of the quadrature rules themselves. To come back to the observation in comment #1: there a half-infinite interval is integrated over. And this is something all these quadrature rules cannot do natively. What they do instead is to transform it (essentially by substitution of the integration variable) to an integral over a finite domain. Of course, there are different possibilities to do that. Probably there is not much logic within the integrators to decide that, they will use just one hard-coded transformation, and possibly the one within quad just by chance gives a very well-behaved transformed integral, while the transformation used by the others doesn't. We are probably not in a position to improve on that, but it would be perhaps nice to state the used transformations in the documentation, so that the user can decide which integrator to use.

Michael Leitner <mleitner>
Thu 05 May 2022 03:07:01 AM UTC, comment #9: 

pushed http://hg.savannah.gnu.org/hgweb/octave/rev/e8ced722b19e that adds the optional output argument to integral and adds a few BISTs for that output.

I'll take a look at the quadcc warning. a quick peek looks like it should be able to check and warn just before the returns at the end.

Nicholas Jankowski <nrjank>
Group Member
Wed 04 May 2022 06:16:04 PM UTC, comment #8: 

I don't see a reason to have a path to calling quad() within integral.m.  The documentation is clear that integral() is a wrapper for various possible quadrature functions.  As long as integral returns a warning when it fails then I think users will understand that they need to go to the individual functions directly.

Rik <rik5>
Group administrator
Wed 04 May 2022 05:49:58 PM UTC, comment #7: 

@jwe, Re: comment #5:  of course. poor wording. 

https://www.mathworks.com/help/matlab/ref/integral.html

shows only

Syntax
q = integral(fun,xmin,xmax)
q = integral(fun,xmin,xmax,Name,Value)


so, since matlab's integral only has one return variable, us adding optional additional sholudn't break compatability (unless they change their interface of course.)

@rik, that sounds like a good plan. so (1) mod quadcc to warn as you said. (2) maybe add warning output option to integral that should get something to pass along from whichever integrator it calls?

is there thought to be any need establish a path by which integral calls quad? would have to make quite a few changes to the desicison tree for that i think, and the recursion issue remains.  integral doesn't have a 'choose your integrator' option as a name/value parameter set, but adding one would be one way to enable quad without playing with selection logic, and it wouldn't break compatability. (i think dblquad and triplequad do that).  At that point though... I'd say just call the one you want directly.

Nicholas Jankowski <nrjank>
Group Member
Wed 04 May 2022 05:33:35 PM UTC, comment #6: 

I would also note that there is a page in the Octave manual about selecting the correct numerical quadrature function based on the characteristics of the input integrand.  'quad' is recommended for Low accuracy and non-smooth integrands.  However, quadcc is general purpose offering Low to High accuracy, smooth or non-smooth integrands, discontinuities, oscillatory functions, finite or infinite limits.  In short, it's a good all-around choice when not much is known about the integrand itself.  Because it is written in C++ it is also fast.

I think we should modify quadcc to check whether it has met its relative and absolute tolerance tests at the end of integration.  If the user has not requested the ERR output then it should issue a warning.  Otherwise, it should stay quiet and assume that since the programmer has requested the ERR output they will be checking for convergence themselves.

Rik <rik5>
Group administrator
Wed 04 May 2022 05:20:39 PM UTC, comment #5: 

By "looking at matlab's integral" I assume you mean looking at the documented interface?

John W. Eaton <jwe>
Group administrator
Wed 04 May 2022 04:52:41 PM UTC, comment #4: 

It's been a while but when I threw integral together I think I looked at the other integral wrappers, dblquad, triplequad, etc., and noticed they used quadcc as their general integrator. I also noticed the restriction mentioned in the help for quad that it cannot be called recursively. since I didn't want to impose that restriction on integral (can recursive calls be identified?), I decided to follow suit.

I was a bit surprised to see that quadcc didn't issue any warning with that poor convergence, noting that quadgk did. I had been under the assumption that the underlying integrater would produce any necessary errors or warnings tied to the actual integration process. Don't know if that itself might be worth looking at.

Also, looking at matlab's integral, there is no compatability issue with adding extra output arguments, so in addition to your suggestion it could make sense to attempt reporting error as an output variable as well. Since the different integrators don't necessarily handle or report error the same way that would be a minor complication, but probalby a surmountable one.

Nicholas Jankowski <nrjank>
Group Member
Wed 04 May 2022 03:38:26 PM UTC, comment #3: 

This is always the case with numerical computing, as opposed to pure math.  You need to check assumptions and verify that the algorithm has worked as expected.  In this case, there are additional outputs from quadcc which give a bound on the absolute error as well as the number of points used for evaluation.  When I call the function requesting this information I get


[q, err, num_pts] = quadcc (f, 1, Inf)
q = -9.7842
err = 2.2304
num_pts = 3101


The absolute error of +/- 2.23 indicates that the algorithm essentially failed since the default absolute tolerance is 1e-10.

I would note that quadgk also has an ERR output.  It might be useful to have integral.m compare this output with output q and issue a warning if they are too divergent.

Another possibility is to have integral also map some integrands to quad.  I don't believe quad can handle complex integrands, but it can handle infinite endpoints.  I did some benchmarking and I don't see much difference between quad and quadcc as they are both compiled functions (the first in Fortran and the second in C++).

Rik <rik5>
Group administrator
Wed 04 May 2022 03:06:29 PM UTC, comment #2: 

fyi rereading the integral function, it only chooses integrators other than quadcc based on parameters being passed that quadcc doesn't support, it doesn't depend on the integral or limits.  so just sending the minimal 3 inputs (function, lowerlim, upperlim) will always call quadcc. specifying waypoints will force it to call quadgk.

Nicholas Jankowski <nrjank>
Group Member
Wed 04 May 2022 02:56:55 PM UTC, comment #1: 

as a comparison, in matlab 2022a:


>> integral(f,1,inf)
Warning: Reached the limit on the maximum number of intervals in use. Approximate bound on error is   2.2e-03. The integral may not exist, or it
may be difficult to approximate numerically to the requested accuracy.
> In integralCalc/iterateScalarValued (line 372)
In integralCalc/vadapt (line 132)
In integralCalc (line 83)
In integral (line 87)

ans =

  -9.993656849721580


So even matlab recognizes some issues with this one.

As the help describes, Octave's integral is a wrapper to quadcc, quadgk, or quadv,selected based on the type of inputs/options chosen, and better accuracy may be obtained by running them directly. quadv breaks on that one, but if I run the other two manually in octave, it appears integral is choosing quadcc:


>> quadcc(f,1,inf)
ans = -9.784172499874394

>> quadgk(f,1,inf)
warning: quadgk: maximum interval count (650) exceeded
warning: called from
    quadgk at line 393 column 7

warning: quadgk: Error tolerance not met.  Estimated error 0.00605827
warning: called from
    quadgk at line 404 column 5

ans = -9.990824621478369

>> quadgk(f,1,inf, "maxintervalcount", 100000)
warning: quadgk: non-finite integrand encountered
warning: called from
    quadgk at line 348 column 7

warning: quadgk: Error tolerance not met.  Estimated error 0.00178496
warning: called from
    quadgk at line 404 column 5

ans = -9.994594924648961


funny enough, the best answer in Octave comes from just calling the 'old' routine 'quad:

Octave:
>> quad(f,1,inf)
ans = -9.999999999999767

matlab:
>> quad(f,1,inf)
Warning: Infinite or Not-a-Number function value encountered.
> In quad (line 100)
ans =
   NaN


it may be that there's some improvement that can be done in function selection, but the goal when implementing integral was to get a compatible function in place. it hasn't really been revisited for improvement since.

Nicholas Jankowski <nrjank>
Group Member
Wed 04 May 2022 02:31:12 PM UTC, original submission:  

The commands

>> format long % Display lots of digits
>> f = @(t) -1./t.^1.1;
>> integral(f,1,inf)


should return 10 but returns 9.784... .  Changing error tolerances does not help. 


Anonymous

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

Attach Files:
   
   
Comment:
   

Attached Files
file #53183:  quadccwarning_bug62412_v2.diff added by nrjank (1KiB - application/octet-stream - adds conditional error warning to quadcc)
file #53177:  quadgk_warn_suppress_bug62412.diff added by nrjank (996B - application/octet-stream - suppresses quadgk final tolerance warning if err return requested)
file #53176:  quadccwarning_bug62412.diff added by nrjank (886B - application/octet-stream)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by mmuetzel (Posted a comment)
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by jwe (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by nrjank (Posted a comment)
  •  

    There are 0 votes so far. Votes easily highlight which items people would like to see resolved in priority, independently of the priority of the item set by tracker managers.

    Only group members can vote.

     

    Follow 13 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-05-15 rik5 StatusReady For Test Fixed
        Open/ClosedOpen Closed
    2022-05-10 rik5 StatusPatch Submitted Ready For Test
    2022-05-09 nrjank Attached File- Added quadccwarning_bug62412_v2.diff, #53183
    2022-05-05 nrjank Attached File- Added quadgk_warn_suppress_bug62412.diff, #53177
    2022-05-05 nrjank Attached File- Added quadccwarning_bug62412.diff, #53176
        StatusConfirmed Patch Submitted
    2022-05-04 rik5 Item GroupNone Missed Error or Warning
        Summaryintegral gives consistent inaccurate result for improper integral to inf quadcc does not warn when quadrature fails to converge
    2022-05-04 nrjank Severity3 - Normal 2 - Minor
        StatusNone Confirmed
        Release6.2.0 dev
        Operating SystemMicrosoft Windows Any

    Back to the top

    Powered by Savane 3.13-4b48.
    Corresponding source code