Test failures due to tolerance in fftfilt.m

[Daniel J Sebald]

I'll toss this one to Ed and Rik, since we were just talking about 
precision issues for svds test failures...

I checked the current state of tests and found this failure:

> > > > > processing /usr/local/src/octave/octave/octave/scripts/signal/fftfilt.m
  ***** test
 r = sqrt (1/2) * (1+i);
 b = b*r;
 assert (fftfilt (b, x  ), r*[1 1 0 0 0 0 0 0 0 0]  , eps);
 assert (fftfilt (b, r*x), r*r*[1 1 0 0 0 0 0 0 0 0], eps);
 assert (fftfilt (b, x.'), r*[1 1 0 0 0 0 0 0 0 0].', eps);
!!!!! test failed
assert (fftfilt (b, r * x),r * r * [1, 1, 0, 0, 0, 0, 0, 0, 0, 0],eps) 
expected
 Columns 1 through 3:
   0.00000 + 1.00000i   0.00000 + 1.00000i   0.00000 + 0.00000i
 Columns 4 through 6:
   0.00000 + 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i
 Columns 7 through 9:
   0.00000 + 0.00000i   0.00000 + 0.00000i   0.00000 + 0.00000i
 Column 10:
   0.00000 + 0.00000i
but got
 Columns 1 through 3:
   0.00000 + 1.00000i   0.00000 + 1.00000i   0.00000 - 0.00000i
 Columns 4 through 6:
  -0.00000 + 0.00000i  -0.00000 + 0.00000i   0.00000 + 0.00000i
 Columns 7 through 9:
  -0.00000 + 0.00000i   0.00000 - 0.00000i  -0.00000 + 0.00000i
 Column 10:
  -0.00000 + 0.00000i
maximum absolute error 2.22478e-16 exceeds tolerance 2.22045e-16
shared variables   scalar structure containing the fields:
    b =
       1   1
    x =
       1   0   0   0   0   0   0   0   0   0
    r = [](0x0)

This is another one of those machine precision issues where I believe 
expectations are too great.  Filtering using the FFT is a numerical 
technique which takes the FFT of both sequences, multiplies the FFTs 
then takes the inverse FFT to get the output sequence.  That's a lot of 
numerical computations, and there must be precision effects.

There are a bunch of tests in fftfilt.m related to tolerance:

%!shared b, x, r
%!test
%! b = [1 1];
%! x = [1, zeros(1,9)];
%! assert (fftfilt (b,  x  ), [1 1 0 0 0 0 0 0 0 0]  , eps);
%! assert (fftfilt (b,  x.'), [1 1 0 0 0 0 0 0 0 0].', eps);
%! assert (fftfilt (b.',x  ), [1 1 0 0 0 0 0 0 0 0]  , eps);
%! assert (fftfilt (b.',x.'), [1 1 0 0 0 0 0 0 0 0].', eps);

The reason the above tests are passing within machine tolerance is that 
they are sort of degenerate cases.  Most of the input values are 0 or 1, 
so the butterflies and accumulations are rather simple.  Those tests are 
fine as is.

It's the next few where precision effects begin to creep in because 
irrational numbers (i.e., number utilizing full precision) are input:

%!test
%! r = sqrt (1/2) * (1+i);
%! b = b*r;
%! assert (fftfilt (b, x  ), r*[1 1 0 0 0 0 0 0 0 0]  , eps);

Again, x = [1 0 0 0 0 0 0 0 0 0], so b convolved with x doesn't entail 
many multiplies except with 0 inside the butterflies.  Again, no issue.

%! assert (fftfilt (b, r*x), r*r*[1 1 0 0 0 0 0 0 0 0], eps);

Here, this one now has some precision issues going on.  There are some 
multiplications in the butterflies probably going one or two levels deep 
(multiplies, sums, etc.) and with numbers utilizing the extent of the 
floating point register.  So, what we are looking at is the effects of a 
couple multiplies.  So, how about a tolerance of 1.5*eps?  2*eps?

%! assert (fftfilt (b, x.'), r*[1 1 0 0 0 0 0 0 0 0].', eps);

This is the same as two previous, just reoriented the x vector--no issue.

ANOTHER ISSUE

Note that the above tests really don't check the integrity of the 
algorithm, because with such degenerate cases the algorithm could easily 
pass even if there were some fundamental flaws.  The test that follows 
the above is meant to address that by using inputs which are all 
non-zero, but the test has a big flaw.  Take a look:

%!test
%! b  = rand (10, 1);
%! x  = rand (10, 1);
%! y0 = filter (b, 1, x);
%! y  = filter (b, 1, x);
%! assert (y, y0);

y0 and y are by definition equal.  I think what is meant is something like:

b  = rand (10, 1);
x  = rand (10, 1);
y0 = filter (b, 1, x);
y  = fftfilt (b, x);
assert (y, y0, 10*eps);

Why pick 10*eps?  Well, I can guess that 5*eps will fail with some 
likelihood.  Between 6 and 10?  I'm not sure.  10 seems like it should 
always pass.  The logic has to be taken with a grain of salt: the 
sequences are length ten, which translates to ceil(log2(10)) = 4 levels 
of butterflies/accumulates.  So that is four levels for which precision 
effects can be introduced.  Possible loss at each level?  Oh, maybe 1.5 
to 2 epsilon on average, given there are a few butterflies per stage? 
So, that is putting us in the range of 6 to 8 epsilon for the overall 
computation.  If one does a bunch of tests with 6*eps:

for i=1:1000
  b  = rand (10, 1);
  x  = rand (10, 1);
  y0 = filter (b, 1, x);
  y  = fftfilt (b, x);
  assert (y', y0', 6*eps);
end

there is failure within a dozen trials.  For 7*eps, there are failures 
within 10 or 20, but also some pretty long runs.  For 8*eps, there is 
failure within 100, but some pretty long runs of 400 or so.  For 9*eps, 
similar to 8*eps.  For 10*eps, failures within 200 trials (I'm already 
proving myself wrong!)...

OK, how about 11*eps?  It's passing 1000 trials most of the time.

So lets go to 12*eps and bump the number of trials up to 1e5...  Got a 
failure!  Way near the end of the run.

We know there is a hard limit on the maximum loss of precision, but 
apparently it is greater than 12*eps.

So, how about 15*eps for a tolerance on that one?  I think that is still 
establishing the integrity of the internal FFT algorithm.  It has passed 
1e6 trials here without failure.  I'm suggesting

%!test
%! b  = rand (10, 1);
%! x  = rand (10, 1);
%! y0 = filter (b, 1, x);
%! y  = fftfilt (b, x);
%! assert (y, y0, 15*eps);

Dan