Optimizing an inner loop calculation in Mathematica

Question

Hi all, I'm currently doing some computation in Mathematica related to Quantum Mechanics. As we've moved from a 1D to 2D lattice model, the problem size is becoming problematic

Currently, we have a summation that looks something like this:

corr[r1_, r2_, i_, j_] = Sum[Cos[f[x1, x2] Angle[i] r1 + f[y1, y2] Angle[j] r2], {x1, HL}, {x2, HL}, {y1, HL + 1, 2 HL}, {y2, HL + 1, 2 HL}];

f[. , .] is a lookup function for a pre-computed correlation function, and Angle[.] is precomputed as well.

There's no way at all to simplify this further in any way. We already took a simple optimization by converting a complex exponential (which had zero imaginary part) to the cosine expression above.

The big problem is that those HL's are based on dimension size: For linear dimension L along an axis, HL corresponds to L^d (d = 2 here). So our computation is O(n^8) in reality, neglecting the sum over i, j.

This normally isn't too bad for L = 8, if it weren't for the fact that we iterate this for 125 values of r1, and 125 of r2 to create an 125 x 125 image.

My question is: How can I most efficiently calculate this in Mathematica? I would do this in another language but there are certain problems that will make it just as slow if I tried it in something like C++.

Extra info: This is a ND-ND (number density) correlation calculation. All of the x's and y's refer to discete points on a discrete 2D grid. The only non-discrete thing here is our r's.

Answer 1

It seems that swapping the Fourier transform with a Cosine transform was the wrong time to optimize, as it hides the fact that this correlation calculation is really just a product of two Fourier transforms (which is the only efficient way to calculate correlations I know of).
With ir1=Angle[i] r1 and ir2=Angle[j] r2 your expression is equivalent to

Sum[Cos[f[x1, x2] ir1 + f[y1, y2] ir2], {x1, HL}, {x2, HL}, {y1, HL+1, 2 HL}, {y2, HL+1, 2 HL}]
== Re@Sum[Exp[I f[x1, x2] ir1] Exp[I f[y1, y2] ir2], {x1, HL}, {x2, HL},{y1, HL+1, 2 HL}, {y2, HL+1, 2 HL}]
== Re[corr1[ir1] corr2[ir2]]

where

corr1[ir_]:=Sum[Exp[I f[x1, x2] ir], {x1, HL}, {x2, HL}];
corr2[ir_]:=Sum[Exp[I f[y1, y2] ir], {y1, HL+1, 2 HL}, {y2, HL+1, 2 HL}];

As I have already cut your scaling exponent in half, I expect you are happy :), but if f is real-valued, you can cut another factor of two of the exponent:
In this case, we can express corr1 as an integral over the values of f -- given that you can somehow get at the weight function w. If nothing else, you can do this numerically with a simple binning procedure.

corr1v2[ir_]:=Sum[ w[fval] Exp[I fval ir], {fval,fvals}],

which makes it clear that corr1 is really just the Fourier transform of the weight function of f (so you should compute it with FFT rather than the sum above). Same goes for corr2.
Alternatively, if f is not real-valued but has enough symmetry to allow you to reparameterize in a form so f only depends on one of the new parameters (say, r,phi), you will also cut down the corr1 integrals to one dimension, although it might not be a simple Fourier transform.