Checking the feasibility of storing plan (in advance) & storing pyfftw.FFTW objects for different input sizes

[NimaSarajpoor]

In PR #25, pyfftw_sdp module was revised. It now stores rfft/irfft object so that it can re-use them later. So, I tried to check a few things:

• How many unique next_fast_n for n in range(2, 2**20 + 1)? Answer: 2623
• How long does it take to perform RFFT planning (with FFTW_MEASURE) for all 2623 input sizes? Answer: 2-3 hours on my machine (macOS M1, 8GB RAM)
• How much is the size of the exported wisdom? Answer: 2.6Mb
• How long does it take to load and import wisdom? Answer: <1 sec
• How much does it take to create all RFFT objects that correspond to those 2623 input sizes? This takes ~3 sec (in total) after importing the wisdom
• How much memory does "2623 stored rfft objects" consume? Answer: TBD

Currently, I am trying to understand the caveats of different memory-profiling approaches. I tried pympler. It shows <1 MB.

[seanlaw]

How long does it take to perform RFFT planning (with FFTW_MEASURE) for all 2623 input sizes? Answer: 2-3 hours on my machine (macOS M1, 8GB RAM)

This sounds like a lot. Should we be concerned? Is this serial or parallel? Can the wisdom independent of architecture?

[NimaSarajpoor]

This sounds like a lot. Should we be concerned?

This is a lot. However, this is a one-time operation. Once they save the wisdom in a file after planning, they should be good to go. Having said that, I understand that users might not be willing to wait that long.

---

Is this serial or parallel?

The planning workload for all 2623 input sizes was performed in a serial manner. As a side, the parameter threads in pyfftw.FFTW class was also set to 1 for all cases.

What if we just play with the number of threads in FFTW, and use less-aggressive flag FFTW_ESTIMATE instead of FFTW_MEASURE? Noticed that MATLAB takes a similar approach by default:

Description
method = fftw('planner') returns the method that the fast Fourier transform functions fft, fft2, fftn, ifft, ifft2, and ifftn use to determine a transform algorithm. The default method is 'estimate', which determines the algorithm based on the size of the data.

---

Can the wisdom independent of architecture?

It is dependent. Section 3.4. in fftw3.pdf says:

There are pitfalls to using wisdom, in that it can negate FFTW’s ability to adapt to
changing hardware and other conditions. For example, it would be perfectly possible to
export wisdom from a program running on one processor and import it into a program
running on another processor. Doing so, however, would mean that the second program
would use plans optimized for the first processor, instead of the one it is running on.

Although it says """it would be perfectly possible to
export wisdom from a program running on one processor and import it into a program
running on another processor"""... my experience was different when I got the wisdom in my machine and tried to use it in Google Colab (will do that experiment again and report back). Maybe OS should be the same but hardware can be different?? Still, according to the document, there might be some performance hit.

[seanlaw]

The planning workload for all 2623 input sizes was performed in a serial manner. As a side, the parameter threads in pyfftw.FFTW class was also set to 1 for all cases.

I was thinking that we can do multiprocessing with single threads and run each wisdom plan independently. It should be embarrassingly parallel. If you have 8 cores, then 2623 might only take 15-20 minutes in total.

What if we just play with the number of threads in FFTW, and use less-aggressive flag FFTW_ESTIMATE instead of FFTW_MEASURE? Noticed that MATLAB takes a similar approach by default:

With the above comment, doing FFTW_ESTIMATE might turn this into under 5 minutes if you use all CPU cores!

[NimaSarajpoor]

The following script uses multiprocessing. The planning of RFFT for all those 2623 cases now takes ~30min. I tested the computed wisdom by loading the file and importing wisdom, and using the flag FFTW_WISDOM_ONLY to make sure the pre-computed wisdoms are used in the process.

import multiprocessing
import pickle
import time

import numpy as np
import pyfftw


def plan_and_get_wisdom(n):
    real_arr = pyfftw.empty_aligned(n, dtype=np.float64)
    complex_arr = pyfftw.empty_aligned(1 + n // 2, dtype=np.complex128)
    pyfftw.FFTW(
        real_arr, complex_arr, direction="FFTW_FORWARD", flags=("FFTW_MEASURE",), threads=1
    )

    # Add IRFFT (the backward transform) here 

    return pyfftw.export_wisdom()

if __name__ == "__main__":
    n_values = [-1]
    for n in range(2, 2**20 + 1):
        next_fast_n = pyfftw.next_fast_len(n)
        # The following if condition may not the most optimal way
        #  but probably better than checking a set for each n
        if next_fast_n != n_values[-1]: 
            n_values.append(next_fast_n)
    n_values = sorted(set(n_values[1:]), reverse=True)
    
    # Create a pool of workers to plan FFTs in parallel
    start = time.time()
    print('Total unique sizes to plan:', len(n_values))
    with multiprocessing.Pool() as pool:
        # wisdom_lst = pool.map(plan_and_get_wisdom, n_values)
        result_iterator = pool.imap(plan_and_get_wisdom, n_values, chunksize=1)
        wisdom_lst = list(result_iterator)
    
    stop = time.time()
    planning_time = stop - start

    start = time.time()
    for wisdom in wisdom_lst:
        pyfftw.import_wisdom(wisdom)
    full_wisdom = pyfftw.export_wisdom()
    with open("wisdom_rfft_only.pkl", "wb") as f:
        pickle.dump(full_wisdom, f)
    stop = time.time()
    merging_and_saving_time = stop - start

    print(f"Planning time: {planning_time:.2f} seconds.")
    print(f"Merging and saving time: {merging_and_saving_time:.2f} seconds.")

[NimaSarajpoor]

How long does it take to perform RFFT planning (with FFTW_MEASURE) for all 2623 input sizes? Answer: 2-3 hours on my machine (macOS M1, 8GB RAM)

This sounds like a lot. Should we be concerned? Is this serial or parallel? Can the wisdom independent of architecture?

If we include the planning of IRFFT in the multiprocessing approach (see my previous comment), the total time will be ~1 hour. So, the total planning time for 2623 cases is reduced from 5-6 hrs to 1 hr.

[seanlaw]

I would do:

num_cores = os.cpu_count()

with multiprocessing.Pool(processes=num_cores) as pool:

because it doesn't make sense to for each process to waste time by switching between the 2K+ wisdom calculations (and possibly using too much memory). Here, I think a process completely finishes one wisdom before moving onto the next.

[seanlaw]

What about using FFTW_ESTIMATE along with multiprocessing? Let's do what Matlab does by default.

At some point, you maybe spending all of your time planning with FFTW_MEASURE than it would take to compute DAMP.

I doubt that many people would wait an hour just to plan when they really want to compute DAMP

[NimaSarajpoor]

I would do:

num_cores = os.cpu_count()
with multiprocessing.Pool(processes=num_cores) as pool:

To be explicit? Want to make sure I am not missing anything else here.

because it doesn't make sense to for each process to waste time by switching between the 2K+ wisdom calculations (and possibly using too much memory). Here, I think a process completely finishes one wisdom before moving onto the next.

Cannot connect the dots... between this part and the code you shared. Can you please elaborate?

Here, I think a process completely finishes one wisdom before moving onto the next.

Right... because, IIUC, each single (worker) process handles the computation for planning for one input size (as I pass one input size to the function), and then, once done, it takes another input size.

---

What about using FFTW_ESTIMATE along with multiprocessing?

With 8 cores, it is less than 15 sec for those 2623 input sizes (including both rfft and irfft). According to fftw3-document.pdf, "FFTW_ESTIMATE" does not do any FFT computation.

"""
The flags argument is usually either FFTW_MEASURE or FFTW_ESTIMATE. FFTW_MEASURE
instructs FFTW to run and measure the execution time of several FFTs in order to find the
best way to compute the transform of size n. This process takes some time (usually a few
seconds), depending on your machine and on the size of the transform. FFTW_ESTIMATE,
on the contrary, does not run any computation and just builds a reasonable plan that is
probably sub-optimal. In short, if your program performs many transforms of the same size
and initialization time is not important, use FFTW_MEASURE; otherwise use the estimate.

You must create the plan before initializing the input, because FFTW_MEASURE overwrites the
in/out arrays. (Technically, FFTW_ESTIMATE does not touch your arrays, but you should
always create plans first just to be sure.)
"""

Let's do what Matlab does by default. At some point, you maybe spending all of your time planning with FFTW_MEASURE than it would take to compute DAMP. I doubt that many people would wait an hour just to plan when they really want to compute DAMP

I agree!

As a side:
Recall that DAMP's algorithm takes advantage of early abandoning, and that basically means it mostly deals with input sizes that are exact power of two. The total planning time for input sizes 2^2 to 2^20 with "FFTW_MEASURE" is less than 20 seconds (with 8 cores). So, we can just do planning for those. And, for any other input sizes that DAMP may encounter, we can go for "FFTW_ESTIMATE". For now, going to park this idea here. We can re-visit this later if needed.

[seanlaw]

Cannot connect the dots... between this part and the code you shared. Can you please elaborate?

Never mind. I incorrectly assumed that pool allowed an unlimited number of "tasks" to be run at the same time. However, it looks like the default processes is the number of cores. Normally, in a Linux kernel, you can have an unlimited number of tasks running within a single process and I was trying to avoid that.

[NimaSarajpoor]

According to our previous discussions, my suggestion for the next steps is as follows:

(1) Change the default planning flag to FFTW_ESTIMATE in pyfftw_sdp.py. Compare the performance of the modified pyfftw_sdp.py with the performance of other modules

(2) Enhance pyfftw_sdp so that it can pre-compute the cutoff input size at which the process should be changed from single-threaded to multi-threaded.

(3) Compare the final version with MATLAB's SDP (Note: In MATLAB, number of physical cores can be controlled via maxNumCompThreads. In pyfftw_sdp, however, the n_threads can be set to the number of logical cores)

[seanlaw]

    n_values = [-1]
    for n in range(2, 2**20 + 1):
        next_fast_n = pyfftw.next_fast_len(n)
        # The following if condition may not the most optimal way
        #  but probably better than checking a set for each n
        if next_fast_n != n_values[-1]: 
            n_values.append(next_fast_n)
    n_values = sorted(set(n_values[1:]), reverse=True)

If you're going to apply set to the list anyways, why not just do:

n_values = []
for n in range(2**20, 1):
    n_values.append(pyfftw.next_fast_len(n))
n_values = list(set(n_values))

This looks much simpler and direct.

[NimaSarajpoor]

@seanlaw
Correct. The way I did it was unnecessarily complicated.

[seanlaw]

Recall that DAMP's algorithm takes advantage of early abandoning, and that basically means it mostly deals with input sizes that are exact power of two. The total planning time for input sizes 2^2 to 2^20 with "FFTW_MEASURE" is less than 20 seconds (with 8 cores). So, we can just do planning for those. And, for any other input sizes that DAMP may encounter, we can go for "FFTW_ESTIMATE". For now, going to park this idea here. We can re-visit this later if needed.

I think this may be a potential opportunity for us to improve the performance of DAMP. DAMP simply chose nextpow2 because it is generally faster. However, for FFTW/pyFFTW, next_fast_len will be faster than simply choosing powers of 2. We should keep this in the back of our minds (i.e., not a priority).

[NimaSarajpoor]

However, for FFTW/pyFFTW, next_fast_len will be faster than simply choosing powers of 2

We can use this idea for the last chunk since original DAMP applies FFT on the Actual number of remaining data points (in the last chunk). Note that we have two different things in DAMP: (1) Checking data points chunk by chunk to take advantage of "early-abondaning" opportunity , and (2) checking last chunk (if no early-abandoning occurs) that contains the remaining data points. I think what you are suggesting here is for part (2). More details below.

---

I feel the following questions/answers were shaped in the head of the authors when they were working on DAMP algorithm:

(1) If I know the best-so-far discord distance "d", and I get a new data point, how should I update that discord distance based on that new subsequence? In other words, how I can figure out if the new subsequence is a discord or not.

(1-1) Just find full distance profile between that new, last subsequence and all historical data. Get maximum value and compare it with "d"

(1-2) Or, let's try some close-to-present data points first (as our first chunk).

Let's go with (1-2). Maybe, in that first chunk, we find a neighbour for that last subsequence and their distance is less than "d". This means we can quickly dismiss the new subsequence as discord and ignore the rest of our search.

(2) Okay. For that first chunk, what should be the size?

2-1) m, the length of subsequence?

2-2) Or, 2m? 10m?

2-3) Or, next_fast_len(m) ?

2-4) Or, nextpow2(m) ?

Let's go with option (2-4) because FFT is generally fast when the input size is exact power of 2.

(3) Okay. Let's say we checked the first chunk and we could not dismiss the last subsequence as discord. What should we do next?

(3-1) Just compute the full distance profile between subsequence and all the remaining data points?

(3-2) Or, just check another chunk, and it's length should be the same as the first chunk.

(3-3) Or, check another chunk but with longer length.

Let's go with (3-3) as it allows us to cover more subsequences in each attempt. If we go with nextpow2(m) in previous step... we can double the length and still the length will be exact power of two.

(4) Okay. So, let's say we check the chunks, one after another, and we have not been able to abandon the search early. Suppose we are now in the last chunk that contains "r" data points. We need to compute the distance profile between that new subsequence (located at the rightmost side of time series) with the remaining data points (located at the leftmost side of time series)

(4-1) Apply fft on array with size r

(4-2) Or, apply fft on array with size nextpow2(r)

(4-3) Or, apply fft on array with size next_fast_len(r)

@seanlaw
Original DAMP goes with (4-1). What you are suggesting is (4-3), which makes sense.

[seanlaw]

@NimaSarajpoor Thank you for providing the extra context (as I am not super familiar with the innards of DAMP)!

We can use this idea for the last chunk since original DAMP applies FFT on the Actual number of remaining data points

Yes, I was thinking about this and realized that it may not be a significant improvement to the performance

Original DAMP goes with (4-1). What you are suggesting is (4-3), which makes sense.

Ultimately, I am thinking about "what is the optimal length of r that will produce the fastest speed before there is a significant increase in computational time?" Sometimes, it might not be next_fast_len(r) but maybe last_fast_len(r) or last_last_fast_len(r). If last_last_fast_len(r) is somehow 5x faster than next_fast_len(r) then it might be more advantageous do last_last_fast_len(r). Do you know what I mean?

I don't recall how r is chosen in the first place. Is it arbitrary?