Computation of triple nested loops as a convolution product?

I'm trying to compute efficiently the following $$$$A_j = \sum_{l'=1}^{\infty}\sum_{k= 0}^{K-1} L_{l'}T_ke^{2\pi i \frac{k}{K}j}\epsilon_{l',k}$$$$ for $$j = 0,1, \ldots, K-2,K-1$$, where $$L_l$$ (dimension $$NL$$) and $$T_k$$ (dimension $$W$$) are 1D vectors of complex numbers and $$\epsilon$$ is a 2D matrix (dimensions : $$L \times W$$) of complex numbers.

Numerically, I'm performing the summation over $$l'$$ from $$1$$ to $$NL$$, where $$N$$ is an integer that has to be chosen sufficiently large (to satisfy conditions unimportant in this context), therefore yielding $$$$A_j = \sum_{l'=1}^{NL}\sum_{k= 0}^{K-1} L_{l'}T_ke^{2\pi i \frac{k}{K}j}\epsilon_{l'\text{mod} L,k},$$$$ a truncated summation over $$l'$$ and the index of $$\epsilon_{l'\text{mod} L,k}$$ being taken with the modulo.

While a straightforward (but naive) implementation of $$A_j$$ in terms of three nested for loops (two for the above double sum and a third that ranges over all $$j$$) works fine, it is terribly slow, especially considering that this implementation is to be repeated for every time step (and there are many of them).

I'm thus looking for an efficient implementation of the above, which looks very much like a 2D separable convolution product. This would allow one to make use of FFT's to compute it quite efficiently, taking advantage that in the frequency domain, a convolution becomes a multiplication of complex numbers.

Does anybody have any idea ?

• What does the $(0)$ attached to $\varepsilon_{l',k}$ imply? Jan 14 '20 at 21:03
• Also, where is the third loop you are speaking of? I see only two. Jan 14 '20 at 21:04
• @WolfgangBangerth, The $(0)$ stands for the initial time and could have been omitted, so that $\epsilon_{l',j}$ are just complex numbers. The third loop comes from the fact that I need to compute $A_j(t)$ for all $j$. Jan 14 '20 at 21:34
• I see. I think that you can make your readers' intuition easier to come if you omit everything that's really not important to your question. This includes the $(0)$ on epsilon, as well as the $(t)$ on all symbols that have such. Simplify the notation until it reflects the core of the problem. Jan 15 '20 at 4:50
• It's somewhat confusing to have L as one of the functions and also the rows of the matrix, and L prime as a summation variable as well doesn't help. But can't you do this first with L (the number of rows of epsilon) FFTs effectively doing the K and J loops, and then summing over the results (the L prime loop) weighing with the L function to get the answer? Jan 16 '20 at 9:10

Let us write it as $$$$A_j = \sum_{l'=1}^{NL} S_{jl'}L_{l'}$$$$

Where

$$$$S_{jl'} = \sum_{k= 0}^{K-1}[T_k\epsilon_{l'\text{mod} R,k}]e^{2\pi i \frac{k}{K}j}$$$$

(I'm using R for the number of rows, the re-use of L is much too confusing)

This last expression is clearly the discrete Fourier Transfrom of the term in the square brackets. Thus evaluting the FT for a given l' gives you the terms for all j in $$S_{jl'}$$. Thus to calculate $$A_j$$ evaluate the FT for all l', and then perform the sum in the first equation, which is simply a matrix vector multiplication.

You can improve this further by exploiting the mod function which means there are in fact only R independent FTs you need to evaluate, as for higher l' the terms will just repeat.

If I have this right this will reudce the number of operations from $$O(NL\times K^2)$$ to $$O(R\times K \times log( K ))$$, assuming the Fourier transform part dominates the matrix vector multiply.