inaccurate results due to inaccurate trigonometric recurrences

[stevengj]

With accurate trigonometric constants, a good FFT algorithm should have $O(\log n)$ error bounds and $O(\sqrt{\log n})$ root-mean-square error, similar to pairwise summation. See also the commentary and references posted with the FFTW accuracy benchmarks.

However, the accuracy of FFTA seems to be much worse, nearly $O(n)$ error growth, presumably because (like many naive textbook presentations) it is constructing the trigonometric constants by a recurrence relation. Here is a plot of the single-precision accuracy on random data (using the double-precision result as the "exact" value for comparison) with both FFTW and FFTA:

via the code:

using LinearAlgebra
relerr(approx, exact) = norm(approx - exact) / norm(exact)

n = 2 .^ (4:20)
err_fft = Float64[]
for n in n
    x = randn(ComplexF32, n)
    push!(err_fft, relerr(fft(x), fft(ComplexF64.(x))))
end
@show err_fft;

Recommendation: the only reasonably efficient way to get accurate trigonometric constants is to precompute them by calling cis (or cispi) and storing the table in the "plan" data structure.

[dannys4]

yeah that's a great point. I'll sit on how to improve this for a little bit---I'd like to allow the user to trade-off between speed/memory/accuracy. My current inclination is to create a choice in the backend using some kind of iterable object and allow you to choose if it computes the values on-the-fly or during the loop. The problem, of course, is type-stability here (which would partially defeat the purpose of keeping both the fastmath-esque code and the more accurate code with precomputation).

In some sense, however, this is a facet of the power-of-two/four case is implemented since it recursively calculates the weights---do you have a sense of how the error would scale if we only precalculate the first twiddle factor? This allows for on-the-fly calculation of the rest of them in-the-loop. It would be a little easier to swallow as a performance/memory bottleneck since we could at least put tuples that scale with $\log(N)$ (as opposed to the memory scaling of $N$ that complete precalculation would require). I also like this because it wouldn't require as much work 😁.

That being said, it's probably best to include a mechanism for scratch space in the nodes, because a similar (but distinct) problem of pre-allocating scratch memory I think came up in #107 .

[wheeheee]

I had also suspected accuracy problems from the twiddle factors when a few tests failed in #107. The w4 = w2 * w2 accidentally made that test pass, but obviously didn't actually solve the accuracy problem.

do you have a sense of how the error would scale if we only precalculate the first twiddle factor? This allows for on-the-fly calculation of the rest of them in-the-loop.

If it's what I think you mean, no. Empirically, doing that is insignificant. I made fft_pow2_radix4! calculate the first twiddle factor with cispi while trying to make tests pass for #107, it did not work. Because most of the error comes from deriving the other twiddle factors via repeated multiplication by w1/w2/w3 in the loop.

Not sure if it's what you want but I think Tasche & Zeuner, 2002 has some other algorithms if you want to make different tradeoffs?

[andreasnoack]

Here is an updated figure after #120