4
$\begingroup$

I'm trying to compute efficiently the following \begin{equation} 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} \end{equation} 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 \begin{equation} 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}, \end{equation} 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 ?

New contributor
HansimGlück is a new contributor to this site. Take care in asking for clarification, commenting, and answering. Check out our Code of Conduct.
$\endgroup$
  • $\begingroup$ What does the $(0)$ attached to $\varepsilon_{l',k}$ imply? $\endgroup$ – Wolfgang Bangerth Jan 14 at 21:03
  • $\begingroup$ Also, where is the third loop you are speaking of? I see only two. $\endgroup$ – Wolfgang Bangerth Jan 14 at 21:04
  • $\begingroup$ @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$. $\endgroup$ – HansimGlück Jan 14 at 21:34
  • $\begingroup$ 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. $\endgroup$ – Wolfgang Bangerth Jan 15 at 4:50
  • 1
    $\begingroup$ 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? $\endgroup$ – Ian Bush Jan 16 at 9:10
2
$\begingroup$

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

Where

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

(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.

$\endgroup$

Your Answer

HansimGlück is a new contributor. Be nice, and check out our Code of Conduct.

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.