# My iterative FFT is slower than my recursive FFT

I've implemented in Python a recursive FFT function and an iterative one. They're both correct (they both pass all unit tests) but, when I plot the complexity curves, I find that the iterative version is slower than the recursive one (the curves are both O(n lg(n)) but the iterative one is shifted up a little). Can you help me identify what, in the following code, is slowing down my iterative FFT method ? Am I using ThreadPoolExecutor wrong ?

def iterative_fft(a: np.array, inverse=False) -> np.array:
"""
:param a: polynomial in coefficient form
:param inverse: direct or inverse DFT flag
:return: DFT(a) in point form
"""
n = len(a)
a = bit_reverse_copy(a)

for s in range(1, int(np.log2(n)) + 1):
m = 2 ** s
omega_m = np.exp(2 * np.pi * 1j / m)
if inverse:
omega_m **= -1

def butterfly(k):
omega = 1
for j in range(m // 2):
t = omega * a[k + j + m // 2]
u = a[k + j]
a[k + j] = u + t
a[k + j + m // 2] = u - t
omega = omega * omega_m

executor.map(butterfly, range(0, n, m))

if inverse:
a = a / n
return a

def bit_reverse_copy(a: np.array) -> np.array:
"""

:param a: polynomial
:return: reordered array based on index bits reverse
"""
n = len(a)
bit_size = int(np.log2(n))
reversed_a = np.zeros(n, dtype=complex)

for i in range(n):
reversed_index = reverse_bits(i, bit_size)
reversed_a[reversed_index] = a[i]

return reversed_a

def reverse_bits(n, bit_size):
"""

:param n: integer
:param bit_size:  bits used for integer representation
:return: reversed bit integer
"""
reverse_n = 0
for i in range(bit_size):
reverse_n = (reverse_n << 1) | (n & 1)
n >>= 1

return reverse_n

• Don't do either? This is an excellent case for never writing this yourself and always using a library. Commented Jul 24 at 16:36

Using threads in a CPU-bound loop is usually a bad idea in Python. Quoting the documentation:

CPython implementation detail: In CPython, due to the Global Interpreter Lock, only one thread can execute Python code at once (even though certain performance-oriented libraries might overcome this limitation). If you want your application to make better use of the computational resources of multi-core machines, you are advised to use multiprocessing or concurrent.futures.ProcessPoolExecutor. However, threading is still an appropriate model if you want to run multiple I/O-bound tasks simultaneously.

The best change you can do to your code to speed it up would be to change the

        with concurrent.futures.ThreadPoolExecutor() as executor:
executor.map(butterfly, range(0, n, m))


call into

        for k in range(0, n, m):
butterfly(k)


or even use a deque from collections to avoid Python's overheads of the for loop:

        deque((butterfly(k) for k in range(0, n, m)), maxlen=0)


This gives a ×5 speedup on my tests.

However, since some indices in range(0, n, m) may get butterfly to go out-of-bounds on arrays whose size is not a power of 2, you may want to shrink your array to a length of 2 ** bit_size or expand it to a length of 2 ** (bit_size + 1) in bit_reverse_copy, if necessary, to ensure no call to butterfly will fail.

Speaking of bit_reverse_copy, making use of numpy vectorized operations can get you a slightly faster computation:

def bit_reverse_copy(a: np.array) -> np.array:
n = len(a)
bit_size = int(np.log2(n))
input = np.copy(a)  # or np.resize(a, 2**bit_size) if you want to shrink it
indexes = np.arange(n)  # or np.arange(2**bit_size) if you want to shrink it
output = np.zeros_like(input, dtype=complex)
reversed_indexes = np.zeros_like(indexes)

for _ in range(bit_size):
reversed_indexes = (reversed_indexes << 1) | (indexes & 1)
indexes >>= 1
output[reversed_indexes] = input
return output


Not seeing the recursive implementation, it is impossible to say why it is faster. That said, few observations:

• ** is a very expensive operation, and should be avoided when possible:

• m = 2 ** s better be m *= 2.
• omega_m **= -1 also better be omega_m = 1.0 / omega_m.
• I would bit a bullet, and implement direct and inverse FFT separately. Testing for if inverse seems like a waste of time.

• Reversed indices may (and should) be precomputed. It doesn't matter in a single invocation, but we usually do FFT more than once. Data vary, indices stay.

It is hard to compare implementations seeing only one. It is possible that you can write the algorithm in a different way that makes it more efficient. Seeing only one specific implementation, all we can do is point out small changes that hopefully reduce cost.

    for s in range(1, int(np.log2(n)) + 1):
omega_m = np.exp(2 * np.pi * 1j / m)
if inverse:
omega_m **= -1


can also be written as

    constant = 2j * np.pi
if inverse:
constant *= -1
for s in range(1, int(np.log2(n)) + 1):
omega_m = np.exp(constant / m)


With this change, you no longer check for inverse inside any loops, so the recommendation in the other answer to write separate forward and inverse functions is no longer relevant.

    for s in range(1, int(np.log2(n)) + 1):
m = 2 ** s


can also be written as

    m = 1
for s in range(1, int(np.log2(n)) + 1):
m *= 2


In butterfly() you write m // 2 three times, twice inside the loop. Compute this once at the start of the function. Same for k + j, which is not an expensive operation, but this is an interpreted language...

        def butterfly(k):
omega = 1
hm = m // 2
for j in range(k, k + hm):
t = omega * a[j + hm]
u = a[j]
a[j] = u + t
a[j + hm] = u - t
omega *= omega_m


concurrent.futures.ThreadPoolExecutor() creates threads inside the loop over s. Creating threads is expensive, consider creating them only once outside the loop.

• Would you double check that np.exp() expression, please, doesn't look quite right. And nit: I think of a constant as a quantity which doesn't change. Consider choosing a different name for the concept of "a quantity whose sign will often toggle". Even a vague factor might work here.
– J_H
Commented Jul 25 at 23:29
• @J_H The sign doesn’t ever toggle. The value is positive for the forward transform and negative for the inverse transform. A negative input to exp is equivalent to elevating the result to the power of -1 (but it comes for free). I called it constant because it is constant inside the loop. I’m sure there’s a better name but I couldn’t come up with it on short notice. Naming things is hard. Commented Jul 26 at 1:37
• That np.exp() edit looks good, thanks. (I started to make that same edit, feeling that I "knew the author's intent", but then the constant name confused me and I waffled on it.) We could make a single assignment to constant using an ... = ... if ... else ... expression, but yeah, that would look clunky. Multiplying by math.copysign(1., inverse - .5) might be a little cleaner.
– J_H
Commented Jul 26 at 1:47
• @J_H Yeah, thanks for pointing that out. I must have copy-pasted the wrong name… :D Commented Jul 26 at 1:49
• @J_H I’m tempted to call the 2π value full_circle. :) Commented Jul 26 at 1:54