bugGNU Octave - Bugs: bug #37613, Octave precision/accuracy is much...

 
 

bug #37613: Octave precision/accuracy is much lower for quadgk

Submitter:  None
Submitted:  Mon 22 Oct 2012 01:59:36 PM UTC
   
 
Category:  Octave Function Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Pramod Kumbhar 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

Wed 25 May 2022 05:57:52 PM UTC, comment #34: 

Michael has already removed this bit of code as part of bug #62468.  I'm going to close this report and the final work can take place there.

Rik <rik5>
Group administrator
Wed 25 May 2022 05:02:16 PM UTC, comment #33: 

Wow, thanks for the diligence rik. It seemingly wouldn't take much to make it size conformant instead of a scalar, but i agree at this point we remove the block. This report will always be here for posterity if a needed use case is later discovered.

Seems the whole function is being reworked over on 62468, should this specific change be made here or in combination over there?

Nicholas Jankowski <nrjank>
Group Member
Wed 25 May 2022 03:59:30 PM UTC, comment #32: 

The discussion is somewhat unwieldly split between this bug report and bug #62468.  I'll comment on the accuracy provision here.

I tried to create a test case that would exercise the bit of code in _quadgk_eval_ that set the "too_close" output.  This turned out to be harder than I thought.  First, I used Michael's corrected expression testing for 100 eps1.  Then, I added a call to warning() so that I would get a sense of how many times this piece of code was activated.


  if (any (abs (diff (x, [], 2) ./ max (abs (x), [], 2)) < 100 * eps1))
  #if (any (abs (diff (x, [], 2) ./ max (abs (x), [], 2))) < 100 * eps1)
    warning ("hit too close");
    too_close = true;
    q = 0;
    err = 0;
    return;
  endif


It turns out, not at all.  I tried with various integrands, and I could not get this code to execute.  Next, I tried cranking up the tolerances to make quadgk work hard to get the correct answer.  I tried setting "RelTol" to 0 which didn't do anything, and so I tried additionally setting "AbsTol" to eps.  At this point the function stopped with a different warning:


warning: quadgk: maximum interval count (650) exceeded


I tried increasing the "MaxIntervalCount" to 1e6, but I was still not able to activate this code.  Finally, I tried directly loosening the criteria for the "too close" condition


  if (any (abs (diff (x, [], 2) ./ max (abs (x), [], 2)) < 1e8 * eps1))


Then with this test case I was finally able to have the code activate.


quadgkv (@(x) sin (x), 0, 1, 'reltol', 0, 'abstol', eps, 'maxintervalcount', 1e6)
warning: hit too close
warning: called from
    quadgkv>__quadgk_eval__ at line 540 column 5
    quadgkv at line 410 column 24

error: q_subs(294656): out of bound 1 (dimensions are 1x1)


But, as you will notice, there was an error.  When the "too close" code is hit it sets "q = 0", but q, in general, should not be a scalar but a size which depends on the number of columns in the "subs" input.

I haven't gone any farther.  This code has never been active in Octave, throws errors when it is activated, and does not have a theoretical basis for existence.  I'm in agreement with Michael that it should be removed.

Rik <rik5>
Group administrator
Thu 19 May 2022 05:07:53 PM UTC, comment #31: 

I was also surprised for a minute now. But the point is that quadgk internally uses a coordinate transform, which transforms also the integrand (multiplying it by the derivative of the old by the new independent variable or something like that), and this broadcasts to the size of the input. Indeed, if you use complex integral boundaries, then no transformation is done, and octave errors out. The next bug found. I will fix it at bug #62468.

Michael Leitner <mleitner>
Thu 19 May 2022 04:43:03 PM UTC, comment #30: 

suprised you were able to run

quadgk (@(x) 1, 1-1e-13, 1)


matlab errors out because the integrators require vectorized functions where size(in) = size(out).  I had to run


quadgk (@(x) ones(size(x)), 1-1e-13, 1)


to test it. Is quadgk not internally vectorized or does octave just do some implicit expansion or something there?

Nicholas Jankowski <nrjank>
Group Member
Thu 19 May 2022 04:06:27 PM UTC, comment #29: 

With respect to my comment #25: at another point in the paper than the one quoted in comment #2, Shampine says "we interrupt the computation when a subinterval is so short that the effects of finite precision arithmetics intrude". Thus, I am very certain that it was never intended that the computation silently proceeds with a contribution to the integral and the error of zero. I would argue to remove the test, as it costs efficiency also when it is false, if it would be true very probably the code will error out with other messages (like a too large number of subintervals), and even if it should finish, the integral will be reasonable.

Michael Leitner <mleitner>
Tue 17 May 2022 08:48:10 PM UTC, comment #28: 

I thought a bit about an interesting example, but didn't succeed, so now a trivial example: if you turn on the test, that is, change the line reproduced below to its intended form

if (any (abs (diff (x, [], 2) ./ max (abs (x), [], 2)) < 100 * eps1))

then

[a,b]=quadgk (@(x) 1, 1-1e-13, 1)

gives an estimate for the integral equal to zero, and also an estimated error equal to zero. If the test is not active, that is, in the present state of the code, you can do the integration for as small intervals as you like -- the estimated value of the integral will be at least in the right order of magnitude (and within a few eps), only the estimated error is strongly over-optimistic. But the latter point is another issue, as already

[a,b]=quadgk (@(x) 1, 1-1e-6, 1)

claims a precision of 7e20, while in reality it is not more accurate than eps(1), which perhaps is a consequence of the transformation.

Michael Leitner <mleitner>
Tue 17 May 2022 03:44:45 PM UTC, comment #27: 

I was trying to construct a case where turning off the line would make a difference (in my expectation, decreasing the accuracy), when I recognized that it is turned off effectively always. But yes, I can now test whether I find a case where turning it on makes a difference. I will do that.

Michael Leitner <mleitner>
Tue 17 May 2022 03:33:06 PM UTC, comment #26: 

comparing the structure of the conditional statement with comment #4, i think you are correct on the parentheses placement error.

Reopening for now so that that can be examined and corrected. would be interesting to see what order of magnitude differences, if any, correcting this line would introduce.

Nicholas Jankowski <nrjank>
Group Member
Tue 17 May 2022 10:18:06 AM UTC, comment #25: 

I am asking to reopen this old report instead of opening a new one, as the relevant discussion is here:

I did some work on quadgk (bug #62468), and while doing that I wondered about the case that sets the too-close condition. First, I do not see why it is necessary to break off at all. If the intervals become so small that you start to see the discretization of the numerical representation of numbers, of course the estimate for the integral of the subinterval will become quite bad (in a relative way). This will stay like that until the subinterval is so small that it covers just two or even one floating-point number, at which point the Gauss and Kronrod estimates will obviously agree, and the subinterval will be accepted (and in my point of view, this is the correct thing to do). More probably, however, is that the max_subintervals limit will be reached, which will give a nice warning, and emit a sensible estimate of the error. So I do not see why Shampine says "We break off the computation ...", and there is no explanation given.

As it is now (seemingly, see below), if the subinterval becomes to small, its contribution as well as its expected contribution to the error is set to zero, and there is no warning. Thus, the algorithm returns as if there was no problem, but it just cuts out those subintervals from the integration range. This will definitely not be what Shampine meant with "We break off the computation...", rather I would interpret this to mean to raise an error.

However: in reality, the code as active since 2016 does not test this. The line in question is


if (any (abs (diff (x, [], 2) ./ max (abs (x), [], 2))) < 100 * eps1)


and if you look closely, you will see that the closing parenthesis of any() is misplaced -- effectively you test whether all left and right entries of the x-table are identical (only then the any () will not return true, that means 0, which indeed is smaller than 100 eps). So actually since 2016 this dropping of subintervals was not applied. There have not been any reports of problems about that, in fact I would expect to have had more reports if subintervals where indeed dropped.

So if nobody has an idea why it is a good idea to "break off the computation" -- and in this case one would have to specify what this should mean, probably not just silently dropping the subinterval -- I would just say to excise this section of the code.

Michael Leitner <mleitner>
Tue 13 Sep 2016 04:05:40 AM UTC, comment #24: 

I finally pushed Marco's patch by way of Lachlan here (http://hg.savannah.gnu.org/hgweb/octave/rev/ed1722f70fad).  Closing the bug from 2012!

Rik <rik5>
Group administrator
Mon 11 Jul 2016 02:43:31 AM UTC, comment #23: 

I forgot to provide an updated patch with the last comment.  Here it is.  (It keeps Marco as the author.)

Marco's is fine without my comments, so take your pick.


(file #37795)

Lachlan Andrew <lachlan>
Mon 11 Jul 2016 02:17:06 AM UTC, comment #22: 

Some trivial feedback on the patch:

  • The new assert in the BIST is much longer than 80 characters.  Could it rewritten to fit under 79 characters?  Perhaps defining



f = @(x) x x .^ 5 .* exp (-x) .* sin (x)


would help.

  • Is it worth retaining the comment "Shampine suggests 100*eps"?


  • The code



+  if (any (abs (diff (x, [], 2) ./ max (abs (x), [], 2))) < 100 * eps)
+    tooclose = true;
+    return;
+  else
+    tooclose = false;
+  endif


hides the fact that this is an early return.  How about


  if (any (abs (diff (x, [], 2) ./ max (abs (x), [], 2))) < 100 * eps)
    tooclose = true;
    return;
  endif

  tooclose = false;


  • Should "tooclose" be "too_close"?


  • I don't like shadowing the function  eps  by a variable.  I know that was in the existing code, but it is also done in the new  _quadgk_eval_.  How about calling it eps1, since it is eps (1)?


Lachlan Andrew <lachlan>
Thu 30 Jun 2016 05:30:37 PM UTC, comment #21: 

Less than 1% change basically means no difference.  That small a time difference could be caused by various background processes waking up on your computer, etc.

Seems like it is worth committing the patch.

Rik <rik5>
Group administrator
Thu 30 Jun 2016 05:21:38 PM UTC, comment #20: 

Actually, I got an average of 0.45156 for the 24 tests of the new version vs. 0.45522 for the 23 tests of the old version.

Marco Caliari <caliari>
Group Member
Thu 30 Jun 2016 03:51:21 PM UTC, comment #19: 

I'm too far out of Octave developpement to be able to give you a definitive yes or no. However, my original objective to this proposed change what the requiring the call to @trans for all of the tests will make it slower. It would be interesting to know how much slower. What does

tic; test guadgk; toc

give before and after your change for exemple ?

David

Anonymous
Thu 30 Jun 2016 02:56:09 PM UTC, comment #18: 

What about the attached patch.

(file #37633)

Marco Caliari <caliari>
Group Member
Thu 30 Jun 2016 02:06:36 PM UTC, comment #17: 

In any case Marco's comments 8 and 10 identified problems in my patch that should be fixed before it is applied

David

Anonymous
Thu 30 Jun 2016 01:47:09 PM UTC, comment #16: 

The issue is not resolved and David's path does not apply anymore. I can try to submit a new patch.

Marco Caliari <caliari>
Group Member
Thu 30 Jun 2016 01:07:32 PM UTC, comment #15: 

Ping, ping, ping!

David, would you like this patch to get in to 4.2.0, or has this issue been resolved?

Lachlan Andrew <lachlan>
Tue 02 Dec 2014 04:06:46 AM UTC, comment #14: 

Ping!

David, do you want to commit the patch you prepared for this bug report?

Rik <rik5>
Group administrator
Mon 17 Mar 2014 09:15:57 PM UTC, comment #13: 

Ping!

David, do you want to commit the patch you prepared for this bug report?

Rik <rik5>
Group administrator
Mon 21 Oct 2013 04:38:23 AM UTC, comment #12: 

David, do you want to go ahead and apply your changeset from last year?

Rik <rik5>
Group administrator
Tue 30 Oct 2012 11:30:01 AM UTC, comment #11: 

I don't think it is "equal to the old one". Please note I changed the variable "x" in _quadgk_eval_ to be "t" and made the new variable "x" be the abscissae in the untransformed space.

Yes it should be the maximum of the absolute value of x.

D.

David Bateman <dbateman>
Group Member
Tue 30 Oct 2012 08:10:19 AM UTC, comment #10: 

I agree with your comment, but your new patch is equal to the old one.
By the way, shouldn't you take the maximum of the absolute value of x?

Marco

Marco Caliari <caliari>
Group Member
Mon 29 Oct 2012 08:48:28 PM UTC, comment #9: 

Opps, yes it should be 'transformed' to the 'untransformed' form with @trans. However, thinking about it, that going to be fairly expensive if we test the distance between all of the abscissae. Also it seems difficult to understand what Shampine meant by


... become closer in a relative sense than 100 times the unit roundoff eps.


and my previous test was wrong for yet another reason; It trying to give a "relative value" I effectively always renormalized the distance between abscissae by the total length between the abscissae (ie 2) so that the value returned in the test would be 0.021 (ie the distance between the 1st and 2nd abscissae divided by 2) so comparing to 100 * eps would always return false.

Thinking about it Shampine probably meant to compare "DeltaX / X < 100 * eps". As the test is relatively artifical we can get the same effective result by comparing between any two abscissae in the interval and so a test like


x = trans ([t(:,1), t(:,end)]);
if (any (abs (diff (x, [], 2) ./ max (x, [], 2))) < 100 * myeps)


Should do the right thing. Here is a modified changset

D.g

(file #26853)

David Bateman <dbateman>
Group Member
Mon 29 Oct 2012 08:28:05 AM UTC, comment #8: 

If 'untransformed' is (0,Inf) and 'transformed' is (0,1), according to Shampine the check should be done in (0,Inf), i.e. with trans(abscissa). So, we should pass @trans to @__quadgk_eval__

Marco

Marco Caliari <caliari>
Group Member
Fri 26 Oct 2012 08:26:55 PM UTC, comment #7: 

Yes I think you're right the simplification will give the same results much faster.

In fact I'm no longer testing the minimum relative size of an interval, but rather the minimum distance between abscissa, since as you said the "arguments" mentioned by Shampine probably mean the distance between abscissa or evaluation points rather than the sub-interval length.

The test is also after the conversion from the transformed interval to the real evaluation points so I'm testing in the interval (a, b).

D.

David Bateman <dbateman>
Group Member
Fri 26 Oct 2012 03:25:49 PM UTC, comment #6: 

Another thing,

are you not testing now in the trasformed interval, i.e. [0,1]?

Marco

Marco Caliari <caliari>
Group Member
Fri 26 Oct 2012 01:08:04 PM UTC, comment #5: 

Couldn't this line


if (any (abs (bsxfun (@(y,z) y ./ z, diff (x, [], 2), (x(:,end)-x(:,1))))(:) < 100 * myeps))


be simplified into


if any(min (abs (diff (x, [], 2)), [], 2) < 100 * myeps * (x(:,end)-x(:,1)))


?

Marco

Marco Caliari <caliari>
Group Member
Thu 25 Oct 2012 09:20:49 PM UTC, comment #4: 

If there is an infinite interval, then h0 is always infinite and the test


if (any (abs (diff (trans (subs), [], 2) / h0) < 100 * myeps))


is always going to suceed and the quadrature exit immediately. I suppose you're right Marco in that matlab must test at the abscissa in the untransformed subintervals rather than do a single overall test on the bounds of the subintervals. Though Shampine doesn't make it clear why the test MUST be done in the untransformed intervals I supposed its safer to do so.

I don't think the changes necessary will be that major. We just need to move the test into _quadgk_eval_ where the abscissa are available. As we would then be dealing with a matrix of abscissa rather than a vector of sub interval end points the test would either need to use bsxfun or the new broadcasting code.

I propose the attached changeset, though a little testing before committing it might be good idea. In any case with this changeset I now see


octave:24> f = @(x)x.^5.*exp(-x).*sin(x);
octave:25> [q,errbnd] = quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12)
q = -15.000
errbnd =  1.0290e-08


So it appears to work correctly.

D.



(file #26829)

David Bateman <dbateman>
Group Member
Thu 25 Oct 2012 08:29:27 AM UTC, comment #3: 

Dear David,

maybe the key is in the word "arguments": you took the length of subintervals as arguments (and in the case [0,inf) there is always a subinterval of infinite length), but in the paper the arguments are the evaluation points for the integrand (see, for instance, the third line in p.4). Since inf is never an evaluation point, the distance between "arguments" is never infinite.
This would require some work on the code...

Marco

Marco Caliari <caliari>
Group Member
Wed 24 Oct 2012 10:05:38 PM UTC, comment #2: 

In fact if you replaced h0 with h you'd also have to replace trans(subs) with subs. As for the sense of the comment, I believe I based it on the statement in section section of the Shampine paper


We break off the computation if arguments become closer in a relative sense
than 100 times the unit roundoff eps. This is not entirely straightforward
because of the transformations: The test must be done after the variable t used
in the formulas is transformed to the variable x given to the function f (x) and
there are four sets of transformations.


The paper itself can be found at

http://faculty.smu.edu/shampine/rev2vadapt.pdf

and I didn't have any other reasons for this comment. However changing line 322 of quadgk to read


if (any (abs (diff (subs, [], 2) / h) < 100 * myeps))


then running "test quadgk" passes all of the tests. These tests include tests for infinite bounds and singular integrands and so probably give a pretty clear indication that such a change would be acceptable, even if Shampine himself didn't recommend it. Perhaps Mathworks modified this behavior after the inclusion of quadgk in Matlab?

D.

David Bateman <dbateman>
Group Member
Wed 24 Oct 2012 12:59:30 PM UTC, comment #1: 

Hi,

instead of removing the check, I think that h0 (length of untransformed interval) should be replaced by h (length of transformed interval). But the comment relative to this check explicitly speaks of "untransformed subintervals". Maybe David had some counterexample in mind.

Marco

Marco Caliari <caliari>
Group Member
Mon 22 Oct 2012 01:59:36 PM UTC, original submission:  

Hello All,

I am new to Octave and running Matlab code on Octave 3.4.3. I am getting incorrect results for my simulations, basically because of lower precision/accuracy of Octave quadrature routine quadgk.

For example, say we have

f = @(x)x.^5.*exp(-x).*sin(x);

In Matlab I get:

K>> [q,errbnd] = quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12)

q      = -1.499999999999836e+01
errbnd =  9.438576028474561e-09

But in the Octave, I get

debug> [q,errbnd] = quadgk(f,0,inf,'RelTol',1e-8,'AbsTol',1e-12)
warning: quadgk: Error tolerance not met. Estimated error 0.0329476
q      = -1.49999941897813e+01
errbnd =  3.29475856720323e-02

From above simple test, I see Octave accuracy is very low.

To find the cause of this, I loooked into the implementation of quadgk routine in Octave (here:http://fossies.org/dox/octave-3.6.3/quadgk_8m_source.html).

The problem is, when one of the limit of integral is Inf, current implementation don't divide subintervals further due to below condition (line no 322):

if (any (abs (diff (trans (subs), [], 2) / h0) < 100 * myeps))
   .....
   break
endif

When a or b is Inf, h0 is Inf and hence code returns results with very low accuracy.

I am not a mathematician and not sure if this is expected behaviour. But simply removing above check we get more accurate results similar to Matlab.

Can someone please help in this issue? Suggestions for getting more accurate results?

Regards,
Pramod.

Anonymous

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #37795:  bug_37613_quadgk.cset added by lachlan (4KiB - text/x-diff)
file #37633:  quadgk.cset added by caliari (3KiB - application/x-tex)
file #26853:  changeset added by dbateman (3KiB - application/octet-stream)
file #26829:  changeset added by dbateman (3KiB - 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 nrjank (Posted a comment)
  • -email is unavailable- added by mleitner (Posted a comment)
  • -email is unavailable- added by lachlan (Posted a comment)
  • -email is unavailable- added by rik5 (Posted a comment)
  • -email is unavailable- added by dbateman (Posted a comment)
  • -email is unavailable- added by caliari (Posted a comment)
  • -email is unavailable- added by None (Submitted the item)
  •  

    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 16 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2022-05-25 rik5 StatusNeed Info Fixed
        Open/ClosedOpen Closed
    2022-05-17 nrjank CategoryLibraries Octave Function
        Assigned tocaliari None
    2022-05-17 nrjank StatusFixed Need Info
        Open/ClosedClosed Open
    2016-09-13 rik5 StatusPatch Submitted Fixed
        Open/ClosedOpen Closed
    2016-07-11 lachlan Attached File- Added bug_37613_quadgk.cset, #37795
    2016-07-11 lachlan Assigned todbateman caliari
        Release3.4.3 dev
    2016-06-30 caliari Attached File- Added quadgk.cset, #37633
    2012-10-29 dbateman Attached File- Added changeset, #26853
    2012-10-25 dbateman StatusNone Patch Submitted
    2012-10-25 dbateman Attached File- Added changeset, #26829
        Assigned toNone dbateman

    Back to the top

    Powered by Savane 3.13-d3ae.
    Corresponding source code