Alternative to (costly) matrix multiplication

Consider the integral

$$2\pi T = -\frac{1}{2} \int_0^{2\pi} A\frac{\partial B}{\partial \xi} \ d\xi$$

where $$A= \sum_{-N}^N i \ sign(n) \ B_n e^{-in\xi}; \quad \quad B= \sum_{-N}^N \ B_n e^{-in\xi}.$$ These expansions imply

$$2\pi T = \sum |k|B_kB_{-k}$$.

Furthermore, $$B_k = \frac{i}{k}\left(-\dot{Y}_k + \sum_{l+j=k}\dot{Y}_lY_j\ j (sign(l)-sign(j))\right),$$

with dots indicating temporal derivatives.

Now, the canonical variables of the system end up being the complex valued variables $Y_{\pm 1},Y_{\pm 2},..,Y_{\pm N}$, and so the Euler-Lagrange equations involve terms like

$$(1) \quad \quad \frac{\partial T}{\partial Y_m}; \quad \quad \frac{d}{dt}\left(\frac{\partial T}{\partial \dot{Y}_m}\right).$$

Therefore, to derive the governing equations, I must explicitly expand terms like these and write them as some function $F(Y_{\pm 1},Y_{\pm 2},..., \dot{Y}_{\pm 1},...)$

Now, the way I'm currently computing this term is to note

$$T = \sum_{k,l} Q_{kl} \dot{Y}_k\dot{Y}_l,$$

where $$Q_{kl} = \frac{1}{4}\sum_{m=1}^N P_{mk}P_{-ml}$$

and $$P_{mn} = (sign(m-n)-sign(n))Y_{m-n} -\frac{1}{2} Y_mY_{-n}.$$

To set up my governing equations via the Euler-Lagrange equations, I need to compute terms like the ones mentioned above in (1). They both lead to the need to calculate

$$\frac{\partial Q_{kl}}{\partial Y_n} = \sum_{m=1}^N \left(\frac{\partial P_{mk}}{\partial Y_n}P_{-ml}+P_{mk}\frac{\partial P_{-ml}}{\partial Y_n}\right).$$

I need to compute this for all n, so I end up needing to perform a series of 2N+1, $((2N+1),N)\times (N,2N+1))$ matrix multiplications. For large $N$, this computation is absolutely slaughtering my run times, and generally take up +80% of the run time for each time step (depending on the size of N). This is my problem.

The way I've been trying to overcome this is to focus my attention on using the best optimized BLAS libraries to compute the matrix multiplications, and exploiting symmetries of the matrices to try to reduce computation time. However, I think I've milked this line of reasoning for all it's worth without getting too deep into the weeds, so I would like to atleast think about alternatives.

I am very new to numerical methods, and might be completely off base, but it seems like there should be a way to manipulate the FFT here, to avoid doing these matrix multiplications as these Cauchy products are similar to convolutions (ie the pseudo spectral method).

However, I don't see how one can explicitly do this, especially as I need to differential these products with respect to the dependent variables of the system.

If anyone has any insight, it's greatly appreciated.

• Your question describes in detail where the problem but omits the key parts that show why exactly you need to do these matrix multiplications :-) Can you rephrase your question by explaining what the problem is you have (not where it comes from) and what you are currently doing to solve it? Commented Feb 25, 2015 at 13:26
• @WolfgangBangerth Thanks for your comment. I have added comments to the body of the text. Hopefully they clarify my point. If not, please let me know. Commented Feb 25, 2015 at 18:04
• How large is $N$ for you? And how many times do you have to do these products? Commented Feb 26, 2015 at 0:11
• It's worth being careful with how you write these equations. For example, if I understand correctly $\partial P_{mk}/\partial Y_k$ is only non-zero for $m=k,2k$, and $\partial P_{-m,l}/\partial k$ is non-zero for $m=k,l+k$. I think your sum has only got four terms in it. The formula for $P$ seems to have wrong indices, $n\leftrightarrow l$, and the index $n$ is missing from the rhs of $\partial Q_{kl}/\partial Y_n$. Commented Feb 26, 2015 at 2:24
• @WolfgangBangerth I have run N up to around 512. A typical example is N=256. Also, for my solutions to remain numerically stable, I need to do this several thousand times. Commented Feb 26, 2015 at 5:36