bugGNU Octave - Bugs: bug #55587, splinefit.m BIST may fail,...

 
 

bug #55587: splinefit.m BIST may fail, slightly out of tolerance result

Submitter:  Mike Miller <mtmiller>
Submitted:  Sun 27 Jan 2019 09:01:08 PM UTC
   
 
Category:  Test Suite Severity:  3 - Normal
Priority:  5 - Normal Item Group:  Inaccurate Result
Status:  Fixed Assigned to:  None
Originator Name:  Open/Closed:  * Closed
Release:  * 5.0.90 Operating System:  * GNU/Linux
Fixed Release:  None Planned Release:  None
* Mandatory Fields

Add a New Comment Rich Markup
   

Jump to the original submission

Thu 31 Jan 2019 04:48:12 AM UTC, comment #11: 

This should be fine for most test runs.

I ran my test script for several hours with the tolerance set to 20⋅eps and still hit a failure case, but it took over 900,000 iterations to find such a case.

Just documenting that this may still fail for users 1 out of a million times, but it's considered good enough for now.

Mike Miller <mtmiller>
Group Member
Thu 31 Jan 2019 12:31:25 AM UTC, comment #10: 

I just ran for an hour or two with this code


more off;

i = 1;
tf = true;

while (tf)
  tf = test ("splinefit");
  i++;
  if (mod (i, 100) == 0)
    printf ("%.3g\n", i);
  endif
endwhile


It got to 314,000 sucessful iterations before I called it quits.

I checked in the change to increase the tolerance to 20 eps here: https://hg.savannah.gnu.org/hgweb/octave/rev/46fa64c0ae40.

Marking as fixed and closing report.

Rik <rik5>
Group administrator
Wed 30 Jan 2019 07:54:55 PM UTC, comment #9: 

The worst tolerance error I have seen is 15.5⋅eps, as shown earlier, after a few hundred or thousand iterations, but I will try increasing to 20⋅eps and run again for a while.

Mike Miller <mtmiller>
Group Member
Wed 30 Jan 2019 07:45:38 PM UTC, comment #8: 

I re-read the original discussion.  I guess we can just relax the tolerance, but I still find it weird that the example you extracted fails repeatedly for you, and yet for Dmitri and I it passes.  That points to something I don't understand which makes me uneasy.

I don't think the tolerance should depend on norm (y).  The test in question calculates expected y values two ways: 1) interp1 and 2) ppval (splinefit ()).  The comparison is done indivdually and independently for each value of y, rather than mixing components of the y vector as would be the case for a solution to Ax = b.

Looking at the code in ppval.m, Octave uses Horner's algorithm to evaluate the cubic polynomial


  for i = 2 : k
    yi .*= dx;
    yi += shiftdim (reshape (Pidx(i,:), dimvec), ndv - 1);
  endfor


In this case, k = 4, so this loop is executed 3 times.  Assuming worst case, every single unit operation results in a discretization error of 1 eps (all of the interpolated values are between -1 and +1 so this is a valid assumption).  That would mean 3 eps for the multiplications and 3 eps for the additions.

After that, you would want to look at the error in the coefficients dx and Pidx. I don't really want to trace down any further in to the algorithms, so we should just pick a number.  Are the dx coefficients off by 2 eps?  Okay that would be 6 eps in total.  Assume the same for the Pidx, and that is another 6 eps.

The assert statement calculates (observed - expected) which could result in another 1 eps error.  I'm at a total of 19 eps so far.

@Mike: Can you change to 20 eps and then see if you ever get any failures?

Rik <rik5>
Group administrator
Tue 29 Jan 2019 03:14:48 AM UTC, comment #7: 

I see that this was last discussed and these tests were changed in 2012.


The test tolerance was increased from 10⋅eps to only 15⋅eps (while also using rand instead of randn). Dan suggested at the time to make the tolerance much larger, between 30⋅eps and 100⋅eps. Ed suggested that the tolerance should be no less than 10⋅eps⋅norm(y), which is about 25⋅eps for the random yb vector.

Mike Miller <mtmiller>
Group Member
Tue 29 Jan 2019 12:56:17 AM UTC, comment #6: 

Mike's specific test case works fine for me (i.e. does not fail)
on i7, Xeon, and amd fx-8350. Is there a way to check if binary is using x87  instead of sse?

Dmitri.
--

Dmitri A. Sergatskov <dasergatskov>
Tue 29 Jan 2019 12:34:40 AM UTC, comment #5: 

Standard build, x86-64 Intel processor, not using --disable-float-truncate.

I've also pasted the 'yb' values from comment #2 into older builds of Octave 3.8, 4.0, 4.2, and 4.4 and they all show the same error for me, "Abs err 3.4417e-15 exceeds tol 3.3307e-15".

Mike Miller <mtmiller>
Group Member
Mon 28 Jan 2019 11:36:21 PM UTC, comment #4: 

Comment #3 should have been "without error".

Rik <rik5>
Group administrator
Mon 28 Jan 2019 10:59:12 PM UTC, comment #3: 

Interesting,  The example in comment #2 runs with error for me on the stable and development branches.

Are you using the option "--disable-float-truncate" during configure?  Are you running an Intel processor?

Also attaching my original test file.

I will use your script to see if I can isolate a failing example of my own.



(file #46116)

Rik <rik5>
Group administrator
Mon 28 Jan 2019 09:18:44 PM UTC, comment #2: 

Here is a specific example that fails for me on both stable and default. I inserted a 'save' into the BIST to get the exact yb vector that causes a failure


>> xb = 0:2:10;
>> yb = [-0.85659524457527947 -0.43282809207909179 -0.82873770262843049 0.85540706303396452 -0.93702483313540752 -0.46131091926622292];
>> x = 0:0.1:10;
>> y = interp1 (xb, yb, x, "spline");
>> assert (ppval (splinefit (x, y, xb, "order", 3), x), y, 15 * eps ());
error: ASSERT errors for:  assert (ppval (splinefit (x, y, xb, "order", 3), x),y,15 * eps ())

  Location  |  Observed  |  Expected  |  Reason
   (101)       -0.46131     -0.46131     Abs err 3.4417e-15 exceeds tol 3.3307e-15 by 1e-16


This was found on the default branch using the following


i = 0;
n = 3;
while (n == 3)
  i++;
  [n, ~] = test ('splinefit.m');
endwhile
i

output:
i =  5361


Mike Miller <mtmiller>
Group Member
Mon 28 Jan 2019 05:46:11 AM UTC, comment #1: 

This works for me on the development branch, but occasionally fails on the stable branch.  I used the following script


N = 1e4;

bm = false (N, 1);

for i = 1:N
  bm(i) = test ("splinefit");
endfor

sum (bm)


On development branch I've repeatedly run this (30,000 attempts) and 0 failures.  On stable I have gotten 1-3 failures in 20,000 attempts.

Rik <rik5>
Group administrator
Sun 27 Jan 2019 09:01:08 PM UTC, original submission:  

I've seen the splinefit.m test fail depending on the random values used in the test:


***** test
 y = interp1 (xb, yb, x, "spline");
 assert (ppval (splinefit (x, y, xb, "order", 3), x), y, 15 * eps ());
!!!!! test failed
ASSERT errors for:  assert (ppval (splinefit (x, y, xb, "order", 3), x),y,15 * eps ())

  Location  |  Observed  |  Expected  |  Reason
   (101)       0.84332      0.84332      Abs err 3.4417e-15 exceeds tol 3.3307e-15 by 1e-16
shared variables   scalar structure containing the fields:

    xb =

        0    2    4    6    8   10

    yb =

       0.76195  -0.60415  -0.95580   0.91181  -0.94831   0.84332

    x =
...


I guess that slightly increasing the tolerance from 15*eps to 20*eps would help.

Mike Miller <mtmiller>
Group Member

 

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

Attach Files:
   
   
Comment:
   

Attached Files
file #46116:  tst_splinefit.m added by rik5 (89B - text/x-matlab)

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by dasergatskov (Posted a comment)
  • -email is unavailable- added by rik5 (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 4 latest changes.

    Date Changed by Updated Field Previous Value => Replaced by
    2019-01-31 rik5 StatusNone Fixed
        Open/ClosedOpen Closed
    2019-01-28 rik5 Attached File- Added tst_splinefit.m, #46116
    2019-01-27 mtmiller Carbon-CopyRemoved 80942 -

    Back to the top

    Powered by Savane 3.13-02a9.
    Corresponding source code