Repository navigation
Speeding up the computation of sliding dot product with fft #938
Description
Activity
@NimaSarajpoor I just noticed that
scipy.fft.rffthas a parameter calledworkers, which can perform FFT in parallel. I wonder if you could try that and see if it makes any difference?Reacted by Nima Sarajpoor@NimaSarajpoor I came across a more recent FFT implementation called OTFFT that claims to be faster than FFTW and has a more generous MIT license. However, I tried to implement the basic
fftfunction in Python but haven't been able to get the same answer asscipy.fft.fft. Here's what I did (List-7: Final version of the Stockham Algorithm):import math import cmath import numpy as np def fft0(n, s, eo, x, y): if not math.log2(n).is_integer(): # Check if n is power of 2 pass m = n // 2 theta0 = 2 * math.pi / n if n == 1: if eo: for q in range(s): y[q] = x[q] else: for p in range(m): wp = complex(math.cos(p*theta0), -math.sin(p*theta0)) for q in range(s): a = complex(x[q + s*(p + 0)]) b = complex(x[q + s*(p + m)]) y[q + s*(2*p + 0)] = a + b y[q + s*(2*p + 1)] = (a - b) * wp fft0(n//2, 2*s, not eo, y, x) def fft(n, x): y = np.empty(n, dtype=complex) fft0(n, 1, False, x, y) for k in range(n): x[k] /= nWould you mind taking a look? Maybe I messed up somewhere but I've been staring at it for too long and I'm not able to spot anything. Thanks in advance!
Reacted by Nima SarajpoorI came across a more recent FFT implementation called OTFFT that claims to be faster than FFTW and has a more generous MIT licene
Cool!
Would you mind taking a look?
Sure! Will take a look.
Also:
I have been tryingscipy.fft.rfft/scipy.fft.fft. Also, as you mentioned before, I am using different number of workers,1vsos.cpu_count(). Haven't seen any improvement yet compared tostumpy.core.sliding_dot_product.According to the scipy doc:
The workers argument specifies the maximum number of parallel jobs to split the FFT computation into. This will execute independent 1-D FFTs within x. So, x must be at least 2-D and the non-transformed axes must be large enough to split into chunks. If x is too small, fewer jobs may be used than requested.
I will test again and share the result and code for our future reference.
Reacted by Sean M. LawAlso: I have been trying
scipy.fft.rfft/scipy.fft.fft. Also, as you mentioned before, I am using different number of workers,1vsos.cpu_count(). Haven't seen any improvement yet compared tostumpy.core.sliding_dot_product.According to the scipy doc:
The workers argument specifies the maximum number of parallel jobs to split the FFT computation into. This will execute independent 1-D FFTs within x. So, x must be at least 2-D and the non-transformed axes must be large enough to split into chunks. If x is too small, fewer jobs may be used than requested.
I will test again and share the result and code for our future reference.
Currently,
core.sliding_dot_productis using thescipy.signal.convolvefunction. In what follows, the performance ofcode.sliding_dot_productis compared with some alternatives.sliding_dot_product_v0 = core.sliding_dot_product def sliding_dot_product_v1(Q, T): n = len(T) X = scipy.fft.rfft(T) * scipy.fft.rfft(np.flipud(Q), n=n) out = scipy.fft.irfft(X, n=n) return out[len(Q) - 1 :] def sliding_dot_product_v2(Q, T): n = len(T) X = scipy.fft.rfft(T, workers=8) * scipy.fft.rfft(np.flipud(Q), n=n, workers=8) out = scipy.fft.irfft(X, n=n, workers=8) return out[len(Q) - 1 :] def sliding_dot_product_v3(Q, T): n = len(T) X = np.fft.rfft(T) * np.fft.rfft(np.flipud(Q), n=n) out = np.fft.irfft(X, n=n) return out[len(Q) - 1 :]And, this is the code for tracking the running time for different window sizes
n = 1_000_000 data = np.array(loadmat('./DAMP_data/mit_long_term_ecg14046.mat')['mit_long_term_ecg_14046'][0]).astype(np.float64) T = data[:n] t = [] for m in range(3, 5000): Q = T[:m] t1 = time.time() comp = sliding_dot_product_function(Q, T) t2 = time.time() t.append(t2 - t1)
As observed:
-
The functions
sliding_dot_product_v1andsliding_dot_product_v2are both usingscipy rfftand they are the same except for the number of workers. As expected, their performances are close to each other. This is because the number of workers affects the performance if we have 2D inputs. -
The performance of
sliding_dot_product_v3(usingnumpy rfft) is close to the existing versionsliding_dot_product_v0.
Reacted by Sean M. Law-
Thanks @NimaSarajpoor. In case it matters (and if you're not already doing this), it would make sense to test window sizes and/or time series lengths in powers of
2rather than increments of1.Reacted by Nima SarajpoorThanks @NimaSarajpoor. In case it matters (and if you're not already doing this), it would make sense to test window sizes and/or time series lengths in powers of
2rather than increments of1.
[Note]
According to the source code ofscipy.fft.rfft, we can use the parameternto pad the inputQwith zeros to make its length the same as the length ofT. so, ifTis power of two, I think we do not need to have power of two for the length of query. I am going to provide the performance of sliding dot product functions for queries with length inrange (4, 1025):
[update] Correction regarding the label of x axis in the bottom figure is: "the length of query "
Reacted by Sean M. Law@NimaSarajpoor I came across a more recent FFT implementation called OTFFT that claims to be faster than FFTW and has a more generous MIT license. However, I tried to implement the basic
fftfunction in Python but haven't been able to get the same answer asscipy.fft.fft. Here's what I did (List-7: Final version of the Stockham Algorithm):import math import cmath import numpy as np def fft0(n, s, eo, x, y): if not math.log2(n).is_integer(): # Check if n is power of 2 pass m = n // 2 theta0 = 2 * math.pi / n if n == 1: if eo: for q in range(s): y[q] = x[q] else: for p in range(m): wp = complex(math.cos(p*theta0), -math.sin(p*theta0)) for q in range(s): a = complex(x[q + s*(p + 0)]) b = complex(x[q + s*(p + m)]) y[q + s*(2*p + 0)] = a + b y[q + s*(2*p + 1)] = (a - b) * wp fft0(n//2, 2*s, not eo, y, x) def fft(n, x): y = np.empty(n, dtype=complex) fft0(n, 1, False, x, y) for k in range(n): x[k] /= nWould you mind taking a look? Maybe I messed up somewhere but I've been staring at it for too long and I'm not able to spot anything. Thanks in advance!
It turns out that
xwill be output if we just avoid dividing it byn.def fft0(n, s, eo, x, y): if not math.log2(n).is_integer(): # Check if n is power of 2 pass m = n // 2 theta0 = 2 * math.pi / n if n == 1: if eo: for q in range(s): y[q] = x[q] else: for p in range(m): wp = complex(math.cos(p*theta0), -math.sin(p*theta0)) for q in range(s): a = complex(x[q + s*(p + 0)]) b = complex(x[q + s*(p + m)]) y[q + s*(2*p + 0)] = a + b y[q + s*(2*p + 1)] = (a - b) * wp fft0(n//2, 2*s, not eo, y, x) # I swapped the params of `fft` function to make its signature similar to `scipy.ftt.ftt` def fft(x, n): y = np.empty(n, dtype=complex) fft0(n, 1, False, x, y) return xAnd, to test it:
for power in range(1, 10): n = 2 ** power x = np.random.rand(n).astype(complex) ref = scipy.fft.fft(x) np.testing.assert_almost_equal(ref, fft(x, n))Reacted by Sean M. LawIt turns out that x will be output if we just avoid dividing it by n.
Hmmm, I wonder why they performed the division?! Thanks for figuring it out. I just ported it over blindly without trying to understand it 🤣.
How about the
ifft?def ifft(n, x): for p in range(n): x[p] = x[p].conjugate() y = np.empty(n, dtype=complex) fft0(n, 1, False, x, y) # for k in range(n): # x[k] = x[k].conjugate()This doesn't seem to match
scipy.fft.iffteither.Would you mind doing a performance comparison if you are able to crack this?
Reacted by Nima SarajpoorIt turns out that x will be output if we just avoid dividing it by n.
Hmmm, I wonder why they performed the division?! Thanks for figuring it out. I just ported it over blindly without trying to understand it 🤣.
How about the
ifft?def ifft(n, x): for p in range(n): x[p] = x[p].conjugate() y = np.empty(n, dtype=complex) fft0(n, 1, False, x, y) # for k in range(n): # x[k] = x[k].conjugate()This doesn't seem to match
scipy.fft.iffteither.Would you mind doing a performance comparison if you are able to crack this?
This should work:
def _ifft(x): n = len(x) # assuming `n` is power of two x[:] = np.conjugate(x) y = np.empty(n, dtype=np.complex128) fft0(n, 1, False, x, y) return np.conjugate(x / n)I am working on some minor changes to boost the performance. I will share the performance of both original version, and the enhanced version, and will compare them with the
core.sliding_dot_product.Reacted by Sean M. LawFor now, I did some enhancements on the new fft / ifft functions suggested in #938 (comment).
Part (I):
I show the performance of four versions against the performance of our reference version, i.e.core.sliding_dot_product. The output of each version is tested to make sure that the function works correclty. The length of time seriesTis$2^{15}$ , and the length of query is inrange(10, 1000 + 10, 10).These are the description of the four versions:
v0 --> use the new functions fft and ifft v1 --> v0 + Reused the already-allocated memory `y` v2 --> v1 + Converted the inner for-loop of `fft` function to a numpy vectorized operation v3 --> v2 + Added njit decorator with `fastmath=True` v4 --> v3 + Parallelized the outer for-loop of `fft` function
For a time series with length
$2^{15}$ , it seems that there is not much difference betweenv0andv1. The changes inv2andv3seem to be very effective. How aboutv4? To better demonstrate its impact, in figure below, I am showing the performance of v3, v4, and the reference only.
And, we can zoom in further by removing the
v3from the figure. Then, we will see:
Part (II):
We now show how thev4performs against the ref (i.e.core.sliding_dot_product) for different length of time seriesT.
As observed, the gap in the performance becomes bigger as the length of time series increases.
The code is available in the notebook pushed to this PR #939.
Next steps:
(1) Make the code cleaner.
(2) Profile the function to see where that increase in the performance gap comes from.
(3) Optimize accordingly.@seanlaw
What do you think?
Also: If I need to add/revise a step, please let me know.Reacted by Sean M. LawIf I understand correctly, the stockham algorithm is NOT faster than
scipy.fft.convolve(or they are about the same after some optimizations). Is that correct? And it also means that the stockham algo is much slower than FFTW?I wonder if there might be some clues in the OTFFT source code in terms of how they might have parallelized it using OpenMP. I'd be pretty happy if we could get within 2x (slower) than FFTW. I'm also guessing that the sawtooth shape observed in
scipy.signal.convolvelikely comes from switching between two FFT algorithms. It's reassuring to see that OTFFT is pretty stable in performance across different distances. I'm confused as to why usingprangewouldn't give you the necessary speedup but that also depends on the hardware that you are using. Maybe I can find some time to test it out on my Mac with many threads.If I understand correctly, the stockham algorithm is NOT faster than
scipy.fft.convolve(or they are about the same after some optimizations). Is that correct? And it also means that the stockham algo is much slower than FFTW?Yes. After some initial optimizations, I can see that the stockham algorithm (from OTFFT) is slower. So, as you mentioned:
- Python Stockham algo (from OTFFT) is slower than
scipy.fft.convolve. scipy.fft.convolveis slower than MATLAB FFTW.
I wonder if there might be some clues in the OTFFT source code in terms of how they might have parallelized it using OpenMP. I'd be pretty happy if we could get within 2x (slower) than FFTW.
I will try to go through it to get some idea. Need to run the tests on MATALB online server if we want to consider MATLAB FFTW as our benchmark.
I'm confused as to why using
prangewouldn't give you the necessary speedup but that also depends on the hardware that you are using.I think it gave us some boost. See the figure below...
v4is the same asv3but with this difference that it usesprange.
numba. get_num_threads()is8in mymacOSsystem.Reacted by Sean M. Law- Python Stockham algo (from OTFFT) is slower than
It appears that maybe we should consider implementing the six step or eight step FFT algorithm next as it should have much better memory locailty and is therefore "optimized". I'd expect this to be faster than our current sliding dot product. I'm not sure how/if any of the radix variants (shown at the same link above) will help.
Reacted by Nima SarajpoorReacted by Nima SarajpoorReacted by Nima Sarajpoorwe should consider implementing the six step or eight step FFT algorithm next as it should have much better memory locailty and is therefore "optimized"
I have implemented six-step-FFT. I will work on eight-step FFT algorithm.
[Note to self]
Before I forget, here are a couple of notes that Imay considerneed to revisit later:- numpy vectorized operation seems to increase the running time (?!)
- Turning off parallelization may speed up the computation (?!)
Reacted by Sean M. LawI have implemented six-step-FFT. I will work on eight-step FFT algorithm.
In case it matters, section 6.3 might be relevant as it discusses Stockham and the 6-step algo. More importantly, it describes how/why cache memory is important
Reacted by Nima Sarajpoor62 remaining items
Also, another thing I have been trying out is to use pre-computed factors instead of computing each factor in fft algorithm. (This is based on what I read before somewhere(?). It is one of the ways one can use to further improve FFT algorithm) The values of these factors ONLY depend on the length of input.
You might be able to use functools.cache (rather than
lru_cache) to memoize the factors:@cache def fft_cosign_factor(i, theta): return math.cos(i * theta) @cache def fft_sine_factor(i, theta): return math.sin(i * theta)So, when you call either of these functions, Python will check whether the given
iandthetacombination has been used before. If not, it will perform the computation and keep the result in memory. If the inputs have been seen before then it will simply pull the result from memory and avoid re-computing the result. Note that this will NEVER work if the input is anumpyarray as it is not possible to ensure that the elements of the array have not been altered. However, in this case, sinceiandthetaare scalar values, then we're good!theta = math.pi / m for k in range(1, m // 2): # factor = math.cos(k * theta) - 1j * math.sin(k * theta) # The next line is equivalent to this factor = fft_cosine_factor(k, theta) - 1j * fft_sine_factor(k, theta)Reacted by Nima Sarajpoor@NimaSarajpoor I discovered something really interesting while looking into cache-oblivious FFT algorithms. The inventor of cache-oblivious algorithms, Matteo Frigo, had provided a cache-oblivious version of FFT in his original paper (see Section 3) and it turns out that Frigo is also the creator of FFTW! Which, you guessed it, uses a cache-oblivious algorithm!
I am starting to think that (in a separate PR) it might be helpful to create a notebook to help explain cache-oblivious algorithms (with diagrams) starting with the basics of how cache sits in between the CPU and main memory (RAM), to understanding blocks, a simple for-loop of how to divide up a matrix and traverse blocks, and how it can be used in a matrix
transpose, and finally to FFT. Ultimately, I'm hoping that this can be generalized to help us think about using blocks/tiles for thestumpfunction as cache-oblivious algorithms (let's call it COA from now on as I'm getting tired of writing it out 😆).Reacted by Nima SarajpoorSo, when you call either of these functions, Python will check whether the given i and theta combination has been used before. If not, it will perform the computation and keep the result in memory. If the inputs have been seen before then it will simply pull the result from memory and avoid re-computing the result.
Precomputing the array has not shown promising result. I am going with the
cacheapproach. Thanks for the very informative explanation!!!Frigo is also the creator of FFTW! Which, you guessed it, uses a cache-oblivious algorithm!
Interesting! I skim the paper before. And, our current cache-oblivious matrix transpose is basically following the approach proposed in the paper. I didn't notice though that he is one of the authors of FFTW! It might be too soon to say this (?) but it seems we are on the right track!
I am starting to think that (in a separate PR) it might be helpful to create a notebook to help explain cache-oblivious algorithms (with diagrams) starting with the basics of how cache sits in between the CPU and main memory (RAM), to understanding blocks,
Agreed! I think that will be fun!
Reacted by Sean M. Law[Update]
I am going with the cache approach.
Apparently the decorators
@cacheand@njitcan work together for a function unless that function itself is being called by another njit-decorated function! (see this GitHub issue)To the see the impact of cache on performance
from numba import njit from functools import cache import math @cache @njit def math_cosine_cache(i, theta): N = 10000000 total = 0 for _ in range(i, i + N): total = total + math.cos(i * theta) return total @njit def math_cosine(i, theta): N = 10000000 total = 0 for _ in range(i, i + N): total = total + math.cos(i + theta) return totalAnd the problem appears here:
@njit def wrapper(i, theta): return math_cosine_cache(i, theta) wrapper(0, math.pi)In fact,
cachecan decorate the njit-decorated function. However,njitcannot decorate cache-decorated function.Reacted by Sean M. LawPrecomputing the array has not shown promising result.
Hmm, this seems surprising. Can you share a small, self contained, reproducible example with only the sine/cosine parts?
Surprising indeed! Will do! Also need to share one to show the impact of different block sizes on performance.
I am trying to see if I can replace the fft recursive function with for-loop. Can that affect the performance? I read that it might be faster as a recursive function pushes each call to stack. I will continue working on it for the next few days and provide an update.
Reacted by Sean M. LawCan you share a small, self contained, reproducible example with only the sine/cosine parts?
I created three version and put them all in this notebook. The first version is the original version. The second version takes advantage of the identity equation
$e^{j\theta} = cos(\theta) + jsin(\theta)$ . In the third version, we pass an array as argument filled with precomputed values. No performance gain is observed in the third version. I think we are dealing with memory access time here.[Update]
Running the code in my MacOS gives me different result. It shows improvement when using precomputed array.Reacted by Sean M. LawAlso need to share one to show the impact of different block sizes on performance.
I checked the impact of blocksize in cache-oblivious algorithm for matrix transpose. You can see this notebook. I do not see any (significant) improvement in changing the blocksize.
Reacted by Sean M. LawI am trying to see if I can replace the fft recursive function with for-loop. Can that affect the performance?
I wrote the for-loop version of the recursive function we use in 6-step FFT (see this notebook). Ran it in MacOS and noticed no performance gain.
Reacted by Sean M. LawI created three version and put them all in this notebook. The first version is the original version. The second version takes advantage of the identity equation ejθ=cos(θ)+jsin(θ). In the third version, we pass an array as argument filled with precomputed values. No performance gain is observed in the third version. I think we are dealing with memory access time here.
Here a few observations for you to consider:
- It looks like
fill_c_arrayisn't even called byfunc_version3so is there even a need to usenjitforfill_c_arrayrather than using straightnumpy? Otherwise, why not considerprangefor filling the contents of the array? - You've split
fill_c_arrayinto two functions (a public function and a private function) and will incur the cost of invoking the private function. This may slow things down ever so slightly (test this with and withoutnjit) but probably not significant enough to matter. Of course, now you are compiling two functions instead of one. - So, when I talking about result caching (i.e., "memoization" with lru_cache), I was really referring to the point that values were constantly being recomputed. However, it looks like the
c_arrayis only being computed once and the values are use only once and that samec_arrayis not being reused for other lengths ofn. This is inefficient. Zooming out a bit and thinking about the bigger picture (DAMP or MADRID), we should remember that we'll be calling FFT many, many times but for different lengths ofn(that are strictly a power of 2), which means that all of our possiblec_arraywould/should havenin this same range too. It makes no sense to precompute the array once and then throw it away and never use it again. Otherwise, we should test for performance of re-using thec_array. Maybe I misunderstood something? - Is the for loop in
_fill_c_arraythe same asnp.power(c_thata, np.arange(m))?
- It looks like
- It looks like
fill_c_arrayisn't even called byfunc_version3so is there even a need to usenjitforfill_c_arrayrather than using straightnumpy? Otherwise, why not considerprangefor filling the contents of the array?
You are right (re: isn't even called by
func_version3)!I initially wrote this function to be used in
func_version3. Later, I decided to not consider it in the running time (as the content of array just depends on its length. It can be computed once and be used in several calls offunc_version3(More on this below)So is there even a need to use njit for fill_c_array rather than using straight numpy? Otherwise, why not consider prange for filling the contents of the array?
I think it is worth it to use
njitto speed it up if our goal is to compute it when DAMP is initialized.You've split fill_c_array into two functions (a public function and a private function) and will incur the cost of invoking the private function.
The private function is actually a recursive function. I can convert it to a for-loop, and then keep everything in one function. I will work on it.
However, it looks like the c_array is only being computed once and the values are use only once and that same c_array is not being reused for other lengths of n.
Two notes:
(1) I simplified the code in the provided notebook. The recursive functionfunc_version3(n, s, x, c_arr)is actually called2 * ntimes by the 6-step FFT algorithm (nis$\sqrt{len(x)}$ ). For instance, for an arrayxwith size2^22, the function is going to be called2 * (2^11)times. So, we are going to use the precomputed array 2^12 times in 6-step function.(2) In DAMP, for example, with each new batch of data point, we need to use FFT. Suppose that the window size
mis set to2^8. When the new batch arrives, we are going to apply FFT on slice with length 2^9, then slice with length 2^10, and so on till we exhaust all historical data point or abandon the search early. Then, we do the same process for the next batch of data. So, in DAMP, for 1M new data points, we may re-use the precomputed array (with length 2^10) 1M times (And, recall that the functionfunc_version3(2^5, s, x, c_arr)is going to be called2 * (2^5)times.)Is the for loop in _fill_c_array the same as np.power(c_thata, np.arange(m))?
Right! Will test its performance as it is definitely cleaner!
Reacted by Sean M. Law- It looks like
I initially wrote this function to be used in func_version3. Later, I decided to not consider it in the running time (as the content of array just depends on its length. It can be computed once and be used in several calls of func_version3 (More on this below)
(1) I simplified the code in the provided notebook. The recursive function func_version3(n, s, x, c_arr) is actually called 2 * n times by the 6-step FFT algorithm (n is ). For instance, for an array x with size 2^22, the function is going to be called 2 * (2^11) times. So, we are going to use the precomputed array 2^12 times in 6-step function.
To properly show the impact of using precomputed array, I can write a function to call
func_version3(n, ...)ntimes (Like what we have in 6-step FFT). This should help us get better / more accurate idea about the impact of this approach on the running time.Reacted by Sean M. Law[WIP]
I initially wrote this function to be used in func_version3. Later, I decided to not consider it in the running time (as the content of array just depends on its length. It can be computed once and be used in several calls of func_version3 (More on this below)
(1) I simplified the code in the provided notebook. The recursive function func_version3(n, s, x, c_arr) is actually called 2 * n times by the 6-step FFT algorithm (n is ). For instance, for an array x with size 2^22, the function is going to be called 2 * (2^11) times. So, we are going to use the precomputed array 2^12 times in 6-step function.
To properly show the impact of using precomputed array, I can write a function to call
func_version3(n, ...)ntimes (Like what we have in 6-step FFT). This should help us get better / more accurate idea about the impact of this approach on the running time.The Colab notebook is now modified to address the point above. I also ran the code in my MacOS and got the following result (code)
In BOTH functions
func_version3andfunc_version4, we are using precomputed array that contains the factors. Infunc_version4, however, we take into account the time required for computing that array.
Now let's consider the full steps 2 & 5 (of our 6-step algorithm) rather than just its cos/sin part (code). I ran it in my MacOS and got the following result:
Note that
fft_v1andfft_v2andfft_v3are exactly the same except for the "factor" part. Infft_v2, we use precomputed array (the approach we took infunc_version4above). Infft_v3, we pass a factor, and, based on that, we compute the other factors (the approach we took infunc_version5above)Reacted by Sean M. LawBtw is
c_thataa typo and should bec_theta?Right!😅 Thanks for catching it! Will fix it.
[Update]
Here, I am providing some important links to external resources or the comments mentioned here:
range(17, 21).Currently, in Stumpy, the sliding dot product [of a query Q and a time series T], is computed via one of the two following functions:
core.sliding_dot_product, which takes advantage of fft trick usingscipy.signal.convolvecore._sliding_dot_product, which uses a njit on top ofnp.dotThe sliding dot product in MATALB (via fft trick) seems to be faster though.
Can we get closer to the performance of MATLAB?