pyFFTW vs MATLAB FFT

[NimaSarajpoor]

pyfftw_sdp is one of the fast methods in computing the sliding-dot-product of two arrays. It uses pyfftw, a wrapper around FFTW, the core algorithm behind MATLAB FFT. So, it is worthwhile to compare pyFFTW and MATLAB FFT to see how close they are to each other from performance perspective.

Note that MATLAB FFT does not support RFFT. Since RFFT is what we care about from sliding-dot-product's perspective, the following comparison is between pyfftw RFFT and MATLAB FFT. The following figures reflect the performance of pyFFTW (RFFT) w.r.t MATLAB (FFT). The timing scripts were executed on MATLAB online, and you can find the scripts attached at the bottom of this comment.

---

The functions, i.e. FFT or RFFT, are called 100 times for each case, and the average of the running time is recorded.

1-Thread pyFFTW (RFFT) vs 1-Thread MATLAB (FFT)

1-, 2-, 4-, 8- Threaded pyFFTW (RFFT) vs 8-Threaded MATLAB (FFT)

It is interesting to see that single-threaded pyFFTW RFFT outperforms MATLAB FFT in almost all cases even if it uses 8 threads.

---

To double check the performance, I increased the number of iterations per case from 100 to 10000. The maximum length of arrays is set to 2^20.

---

Scripts used for timing:

run_FFT_MATLAB.zip
run_pyfftw_on_MATLAB.zip

[seanlaw]

It is interesting to see that single-threaded pyFFTW RFFT outperforms MATLAB FFT in almost all cases even if it uses 8 threads.

Something tells me that using multiple threads is only advantageous when we have more dimensions in our input array.

So, @NimaSarajpoor , what conclusions should be drawn from these results?

[NimaSarajpoor]

Something tells me that using multiple threads is only advantageous when we have more dimensions in our input array.

Or, when the length of array is large. For instance, see the red plot in my previous comment. But, even in that case, single-threaded pyFFTW seems to be performant.

So, @NimaSarajpoor , what conclusions should be drawn from these results?

From my perspective, it confirms that the pyFFTW block performs well. I've started comparing MATLAB-FFT-based sliding-dot-product and pyFFTW-based sliding-dot-product on MATLAB Online's server. If that block shows a good performance too, then my next step is to go for the DAMP block.

[seanlaw]

Note that MATLAB FFT does not support RFFT. Since RFFT is what we care about from sliding-dot-product's perspective, the following comparison is between pyfftw RFFT and MATLAB FFT

Doesn't this mean that MATLAB FFT should be around twice as slow as pyfftw RFFT (I recall some trick that you applied before)? Is this the right comparison? I guess, I'm fine with it as long as the MATLAB+DAMP code also used MATLAB FFT and not RFFT

[NimaSarajpoor]

Doesn't this mean that MATLAB FFT should be around twice as slow as pyfftw RFFT

Asymptotically speaking, yes. Not sure what input size is considered as "large" but if we look at this image from my earlier comment (the first figure), we can see that pyFFTW (RFFT) is 2x/3x faster than MATLAB FFT when the length of array is 2^20 or more (Why 3x?? Probably we do not need to think about it for now. See below)

Is this the right comparison? I guess, I'm fine with it as long as the MATLAB+DAMP code also used MATLAB FFT and not RFFT

Damp_2_0.m (see the attached file, obtained from DAMP's supporting website), uses fft to compute sliding dot product. So, as long as pyFFTW (RFFT) is faster than MATLAB FFT, then we should be good. For now, we can focus on the next block, i.e. sliding-dot-product.

DAMP_2_0.m.zip

[NimaSarajpoor]

So, @NimaSarajpoor , what conclusions should be drawn from these results?

From my perspective, it confirms that the pyFFTW block performs well

Incorrect.

My assumption was that: "Since pyFFTW's RFFT outperforms MATLAB's FFT, we should get a similar conclusion for their corresponding inverse functions." However, that is not correct. Although pyFFTW's RFFT outperforms MATLAB's FFT, this is NOT true for the inverse of those functions. According to the following result, MATLAB's ifft outperforms pyFFTW's irfft when input array is small. In the image below, the x-axis shows the log2 of len(T), where T is the output of pyFFTW's rfft or MATLAB's ifft. The timing of each case is the average of 10_000 iterations.

---

Note: The MATLAB and python scripts will be attached to this comment.

[seanlaw]

According to the following result, MATLAB's ifft outperforms pyFFTW's irfft when input array is small.

Wow, that is interesting and I wonder how that is possible?! It is more than twice as fast in some cases.

[NimaSarajpoor]

It is interesting indeed. The strange performance of pyfftw-based sliding-dot-product got my attention, and encouraged me to check my assumptions.

I do not know at this point. We can see that IRFFT is faster than IFFT for longer arrays. That means IRFFT does not just apply IFFT on full-length array. I may need to check the code flow in pyfftw and see what blocks of codes are triggered.

[NimaSarajpoor]

For complex arrays with no hermitian symmetry, pyFFTW's ifft outperforms MATLAB's ifft. In our case, however, IFFT is applied to array with that symmetry, and it turns out MATLAB's ifft performs considerably better in that case for arrays with lengths <= 2 ^ 8.

In the figure below, we can see that MATLAB’s IFFT performs about 4×-6x faster on a complex array with Hermitian symmetry compared to one without it. (I do not see such optimization in pyFFTW's irfft for arrays with length in that range.)

---

[Update]
Created an issue on pyFFTW's GitHub repo: pyFFTW/pyFFTW#425

[seanlaw]

In our case, however, IFFT is applied to array with that symmetry

Is it because our array is "real"?

[NimaSarajpoor]

@seanlaw

Is it because our array is "real"?

That is correct. From sliding-dot-product standpoint, we have two real-valued arrays T and Q.

• MATLAB: T, Q --> F1 = FFT(T), F2=FFT(Q) --> IFFT(F1*F2)

• pyFFTW: T, Q --> R1 = RFFT(T), R2=FFT(Q) --> IRFFT(R1*R2)

In MATLAB's case, F1 and F2, each, has the hermitian symmetry. Therefore, F1*F2, the input array for MATLAB's IFFT, will have that symmetry too.

[NimaSarajpoor]

[IMPORTANT UPDATE]
FYI: @seanlaw
It seems that the reason for the performance boost shared below is that the thin wrapper .execute() does NOT normalize the output. So, basically, its output is different! Read more here.

As a side note: If we want to just select a certain set of values from the sliding-dot-product (e.g. top-k smallest ones), we can first use the un-normalized output to select those elements and then scale those elements only. IIRC in DAMP, however, we need to update a matrix profile by calculating distance profiles with different lengths, meaning we need to consider that 'normalization' part. To not confuse ourselves, let's try and see if we can have that normalization and get performance boost from other avenues.

---

As suggested by @hgomersall, calling .execute(), the thinnest wrapper around actual execution, can make sure there's no inadvertent copies etc. The result confirms this and we can see a significant boost!! For our convenience , going to share the result below:

You can use this code to obtain the result above. For the thinwrapper case, just replace fftw_obj() with fftw_obj.execute().

So:

import pyfftw

# R is the output of `RFFT(T)`

a = pyfftw.empty_aligned(len(R), dtype=np.complex128)
fftw_obj = pyfftw.builders.irfft(
        a, overwrite_input=True, n=len(T), threads=1
)

a[:] = R
fftw_obj.execute()

[NimaSarajpoor]

Going to close this as this topic has been explored and discussed in #27.