I need to calculate the following linear convolution

$$ y[n] = h_1[n] \circledast h_2[n] \circledast h_3[n] \circledast \cdots \circledast h_{k-1}[n] \circledast h_k[n]$$

where $k$ exceeds $5000$ and the length of each $h_i[n]$ exceeds $100,000$. I tried the method which is based on the circular convolution/DFT but it does not work because each $h_i[n]$ must be zero padded so DFTs then takes too much time and memory. Are there any other algorithm available?

  • 2
    $\begingroup$ what is the meaning of "$n$ exceeds 5000" when the length of $h_i[n]$ already exceeds 100,000? why not just say $n$ exceeds 100,000? $\endgroup$ – robert bristow-johnson Aug 3 '16 at 11:52
  • 1
    $\begingroup$ I need to convolve 5000 sequences, each consistsis of 100,000 elements. I edited question $\endgroup$ – zer0hedge Aug 3 '16 at 12:00
  • 1
    $\begingroup$ @zer0hedge: When implementing via the DFT, are you trying to do one really long zero-padded DFT of each sequence? It's not surprising that you would have problems with that, as the numbers you gave suggest you would need ~500-million-point DFTs. Instead, fast convolution techniques like overlap-save should be applicable to your problem; with this technique, the DFT size you use is only dependent on the lengths of the individual $h_i[n]$ signals, not the intermediate convolution results that grow in length after each step. $\endgroup$ – Jason R Aug 3 '16 at 13:08
  • 1
    $\begingroup$ if you were to do the convolution in the most straight-forward naive way, how much computation would that take? are all of the $h_i[n]$ the same length? $\endgroup$ – robert bristow-johnson Aug 4 '16 at 2:29
  • 1
    $\begingroup$ @MarcusMüller Due to sparsity and the large number of multiplications involved, a leading and trailing parts of the 'running convolution' quickly become zeros. So I do a 'naive' convolution (i.e. apply $h_i[n]$ sequentially), but in addition I cut tails and heads of intermediate results and memorize how many zeros I cut from the head in order to add them in the end to the final result. It gives significant performance improvement, especially during the first steps. It currently takes less than 1 hour to calculate without parallelization. $\endgroup$ – zer0hedge Aug 4 '16 at 14:12

So, I'm becoming a bit notorious with the whole "by using commonly available tools, how complex is your problem really?" business, but oh well:

I went ahead and simulated the "last" of your convolutions, i.e. I went ahead and picked up a FFT FIR (which is convolution by the overlap-save method Jason R mentioned) by dumping 500 Million random samples through a filter of 100,000 taps:

Flow graph

So, this flow graph produces an average total of 42 MS/s, or 420 (500x100,000) convolutions per second. The complete program, which includes initial FFT of the taps, calculating of the taps and the random pool, loading of libraries and set up of all, takes about 12 seconds. So, yes, doing this exact operation 5000 times would exceed the average coffee break's duration, and take nearly 17h.

I tried timing the initial FFT of the tap vector (which the FFT filter blocks automatically), but I failed - it took much shorter than producing the taps vector itself. So I made a separate flow graph that just does 100,000-transforms of input data, and includes all the data generation/copying-overhead in my speed estimate. With that flow graph, my PC does 1900 FFTs per second, so it should take roughly half a millisecond to do a single one, or 2.5 s to do 5,000 of those.

Which also means that if you're smart, and just used the overlap-add algorithm as contained in the FFT FIR block and hoped that performance scales roughly linear with length of the input vector you're transforming, then every convolution would take about $\frac{12}5\text{ ms}= 2.4\text{ ms}$ longer than the previous, giving us a rough estimate of $\sum\limits_{n=0}^{5000} n\cdot 2.4\text{ ms}\approx 8\text{ h},\,25\text{ min}$.

Obviously, that leaves a bit of headroom for optimization. Now, Operlap-Add and Overlap-Save should actually scale exceptionally well to multi-core/distributed solutions. Also, doing a set of long FFTs also sounds like something that GPUs are great at.

| improve this answer | |
  • 1
    $\begingroup$ A direct performance gain is achieved by dividing those 5000 or so cascade convolutions into parallel pairs. Such as $y[n] = h_1[n]*h_2[n]*h_3[n]*h_4[h]$ , brutally requires 3 serial convolutions. But if you first perform $y_1[n]=h_1[n]*h_2[n]$ and $y_2[n]=h_3[n]*h_3[n]$ and finally $y[n]=y_1[n]*y_2[n]$ . This requires 2 parallel and 2 serial convolutions. Similarly 512 such signals would require at most $\log_2(512) = 9$ serial steps instead of 511... of course assuming you can parallely convolve 256 pairs, which is only possible with GPU, a CPU would at best do 8 or 16 ? $\endgroup$ – Fat32 Aug 4 '16 at 12:15
  • $\begingroup$ @Fat32 pretty much what I've been thinking about; a "binary convolution tree". Problem really is that these kind of things vastly benefit from getting things from memory caches,so wildly jumping through "distant" payloads might offset the theoretical reduction in complexity at some point. Also,you're right,a CPU couldn't do so much in parallel (we're often also bound by mem bandwidth), and some of the operations in a single convolution already scale well on multiple CPU cores (FFTw has multi-threading, and point-wise multiplication scales excellently)->parallelize within conv, not across convs $\endgroup$ – Marcus Müller Aug 4 '16 at 12:31

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.