[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Test failures due to tolerance in fftfilt.m
From: |
Daniel J Sebald |
Subject: |
Test failures due to tolerance in fftfilt.m |
Date: |
Thu, 06 Sep 2012 14:59:33 -0500 |
User-agent: |
Mozilla/5.0 (X11; U; Linux x86_64; en-US; rv:1.9.2.24) Gecko/20111108 Fedora/3.1.16-1.fc14 Thunderbird/3.1.16 |
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
- Test failures due to tolerance in fftfilt.m,
Daniel J Sebald <=