Repository navigation
What is the optimal size of block in the overlap-add method? #37
Description
Activity
I know you are referring to the tutorial but, for the actual code,Is it possible/sufficient to simply call a SciPy function to obtain the optimal block size?
Is it possible/sufficient to simply call a SciPy function to obtain the optimal block size?
It is, and this issue should be treated as low-priority.
That scipy's function is _calc_oa_lens and its logic is actually used in PR #35 (as a side note: better to just use that function directly if possible rather than using its pieces. Need to check impact.)
The scipy's approach is based on the Eq.3 in the overlap-add method page on wikipedia. And, analyzing that "equation" might unlock an opportunity for boosting the performance.
Regarding the size of block in overlap-add:
That scipy's function is _calc_oa_lens and its logic is actually used in PR #35 (as a side note: better to just use that function directly if possible rather than using its pieces. Need to check impact.)
It turns out that using the existing function _calc_oa_lens results in a considerable performance hit. The reason seems to be related to the fact that this function does not consider the type of input (real or complex) when computing
next_fast_len. The impact of passing the defaultreal=Falseinblock_size = sp_fft.next_fast_len(math.ceil(opt_size))is shown in the figure shared in this commentReacted by Sean M. LawThe Scipy's function _calc_oa_lens calculates the optimal size of block,
N, in overlap-add method and it is based on the closed-form solution that minimizes the following objective function, which is the computational load per output sample. See Eq3 in wikipedia$$\frac{N \times (log2(N) + 1)}{N - m + 1}$$ I think it is worth it to understand where this equation comes from, and whether it should be revised or not.
Now, let's revisit the problem: we would like to compute the valid convolution between
QandTvia overlap-add method. Suppose we createnbblocks, each with sizeNand supposeNis a power of two. Let's provide an overview of the required computations for overlap-add method using thosenbblocks.(0) Load/Reshape T into the
nbblocks, each with sizeN
(1) Compute FFT of zero-padded Q, where the final length isN-->1time
(2) Compute FFT for a block of sizeN-->nbtimes
(3) Multiply FFT computed in (1) to FFT of each block computed in (2) -->nbtimes
(4) Compute IFFT after multiplication for each block -->nbtimes
(5) Reshape arrays and slice to extract the valid convolutionBefore we talk about the computational load of those steps, let's pay attention to a few notes:
(i) According to FFT page in wikipedia:In contrast, the radix-2 Cooley–Tukey algorithm, for N a power of 2, can compute the same result with only
$\frac{N}{2}log2(N)$ complex multiplications (again, ignoring simplifications of multiplications by 1 and similar) and$Nlog2(N)$ complex additions(ii) According to overlap-add method, the following statement is provided for
N, the size of the block:... such that
N = ..is an integer power-of-2(iii) Also, According to the description of Eq. 3 in overlap-add method:
When the DFT and IDFT are implemented by the FFT algorithm, the pseudocode above requires about N (log2(N) + 1) complex multiplications for the FFT, product of arrays, and IFFT.[B] Each iteration produces
N-m+1output samples, so the number of complex multiplications per output sample is about: .... (Eq. 3)Therefore, it seems that
Nis power of two, and ONLY the complex multiplications are considered for the computational load when deriving Eq 3. Let's check the computational load for the step 0-5 provided earlier by only considering the "complex multiplication":(0) negligible w.r.t to FFT's load
(1)$\frac{N}{2}log2(N)$
(2)$nb \times \frac{N}{2}log2(N)$
(3)$nb \times N$
(4)$nb \times \frac{N}{2}log2(N)$
(5) negligible w.r.t to FFT's loadThe total computational load
Cis:$$C = 0 + \frac{N}{2}log2(N) + nb \times \frac{N}{2}log2(N) + nb \times N + nb \times \frac{N}{2}log2(N) + 0$$ $$C = \frac{N}{2}log2(N) + nb \times N(log2(N) + 1)$$ Recall that each block (of size
N) is padded withm - 1zeros, wherem=len(Q). Therefore, each block containsN - (m - 1)data points fromT. Therefore the number of blocks,nb, is:$$nb = \left \lceil{\frac{n}{N-m+1}}\right \rceil$$ where
n = len(T).Hence:
$$C = \frac{N}{2}log2(N) + \left \lceil{\frac{n}{N-m+1}}\right \rceil \times N(log2(N) + 1)$$ Note that the first term, i.e.
$\frac{N}{2}log2(N)$ , is related the one-time FFT computation of the zero-padded Q. If we ignore the first term, then we almost get (Eq. 3). If we ignore the ceiling operator, then we can dropnsince it is fixed and has no impact on the optimalNthat minimizes the objective function provided in (Eq. 3).Also, what is the domain of the block size
N? It should be in range[m, next_fast_n)The minimum ismbecause it needs to have at least one data point fromTfollowed bym-1zeros (see Tutorial_FFT_Based_SDP.ipynb).Nshould be<next_fast_n. Otherwise, there should be no benefit in splittingTinto blocks since the fft-based circular convolution (direct approach) applies fft on two arrays, each with sizenext_fast_len(n).We are now at a good position to think about a few questions:
- Should we ignore the first term? Can we still get a closed-form solution if the first term is included?
- If we have real-valued arrays, we can use RFFT/IRFFT instead of FFT/IFFT respectively. Can we revise the total computational cost
Caccordingly? - While this is a nice way of representing the total computational cost, we should acknowledge that there exist a couple of caveats: (a) The total cost is only based on complex multiplication, which may not be the dominant factor when
Nis small. (b) The cost is based on Cooley-turkey radix-2 approach, which is forNthat is exact power of two. However, Scipy doesnext_fast_lenonN, which may result in sub-optimal answer.
Reacted by Sean M. LawIn my previous comment, I've obtained a cost function that (approximately) shows the computational cost of overlap-add via circular convolution. This is to provide an overview of the relevant information regarding the computation of valid convolution between Q and T.
m: length of Q n: length of T N: size of block (in overlap-add)(1) cost function of fft-based circular convolution (direct fft approach without overlap-add) :
$3\frac{n}{2}log2(n) + n$
(2) cost function of overlap-add acc to wikipedia:$C = \frac{N(log2(N) + 1)}{N - m + 1}$
(3) optimal length acc to Scipy: closed-form solution of (2), followed bynext_fast_len
(4) optimal length acc to Scipy (enhanced): closed-form solution of (2), followed bynext_fast_len(..., real=True)
(5) cost function of overlap-add acc to my previous comment:$C = \frac{N}{2}log2(N) + \left \lceil{\frac{n}{N-m+1}}\right \rceil \times N(log2(N) + 1)$ If we consider real-valued arrays, RFFT/IRFFT will be used. So, for array with size
x, their computational cost should be close to the computational cost of FFT/IFFT on array with lengthx/2. Note that there will be no change in the number of blocks or elements in the final output array. Therefore:(6) cost function of overlap-add acc to wikipedia with RFFT/IRFFT:
$C = \frac{N/2(log2(N/2) + 1)}{N - m + 1}$
(7) Optimal length acc to Scipy (enhanced) for cost function (6): TBD?
(8) cost function of overlap-add acc to my previous comment but for RFFT/IRFFT:
$$C = \frac{N/2}{2}log2(N/2) + \left \lceil{\frac{n}{N-m+1}}\right \rceil \times N/2(log2(N/2) + 1)$$ (9) Optimal length for the cost function in (8): TBD?
Okay... I think we are probably heading in the right direction! I checked the plots in #35 (comment) and took a closer look at a couple of cases where the
oaconvolvemethod was outperformed by a direct, FFT-based convolution (pocketfft_r2c_c2r_sdp.py).For instance, let's look at the case where
m = len(Q) == 2**16andn = len(T) == 2**20, whereoaconvolveis outperformed. This suggests something might be wrong with its block size!- The objective function is
$\frac{N (\log_2(N) + 1)}{N - m + 1}$ (see Eq. 3 on Wikipedia), and the value that minimizes this is1_017_497. - SciPy's
oaconvolvecomputes the block size by passing this optimal value tonext_fast_len(..., real=False), which yields1_018_710. - In PR Fixed #34 Customizing scipy's oaconvolve #35, however, we pass the value to
next_fast_len(..., real=True), which yields1_024_000. -
[NEW] Let's check out the computational cost (note that we are using RFFT/IRFFT):
-
No overlap-add:
$3 \times \frac{n/2}{2} \times \log_2(n/2) + n/2$ . Forn == 2**20, the cost is15_466_496. -
With overlap-add:
$\frac{N/2}{2} \log_2(N/2) + \left\lceil \frac{n}{N-m+1} \right\rceil \times \frac{N}{2} (\log_2(N/2) + 1)$ . Forn == 2**20,m == 2**16, andN == 1_024_000, the cost is25_300_204, which is higher than the cost of the direct approach shown in (4-1).
-
No overlap-add:
The approximate cost function used in (1), i.e.,
$\frac{N (\log_2(N) + 1)}{N - m + 1}$ , assumes standard FFT/IFFT are used. For RFFT/IRFFT, this cost function becomes$\frac{(N/2) (\log_2(N/2) + 1)}{N - m + 1}$ .Although this cost function is an approximation (a more accurate one was provided in 4-2 above), it has a closed-form solution. We can compute the value that minimizes it, and then pass that value to
next_fast_len(..., real=True)to find the convolution block size. Ultimately, we should plug these values back into (4-1) and (4-2) to verify whether the computational cost of overlap-add is actually lower than the direct approach.Reacted by Sean M. Law- The objective function is
In Tutorial_FFT_Based_SDP.ipynb. the overlap-add method is explained. This method performs several independent circular convolution on the blocks of array. However, it does not provide any guidance on how to find an optimal size of block. It should be a good idea to do an analysis.