# Speed up Mathematica code involving Convolve

I'm trying to plot the averaged autocorrelation function of a random wire profile $\xi(x)$ which is build from unitsteps which are each L units long. To calculate the average I'm generating nMax realizations of the profile, then I proceed to calculate the autocorrelator $\langle \int_{-\infty}^{\infty} \xi(x)\xi(x-y)~dx \rangle$. Brackets indicate averaging over the nMax copies of the profile. My final goal is to plot this expression as a function of y.

I started out by using both the built-in Mean and Convolve functions, however it is ridiculously slow on my machine (late 2007 MacBook). To give you an idea, it takes 150 seconds to calculate corr[0] alone (even with nMax=1).

Is there anything I'm completely missing here? Or can someone of you even hint at a faster and more elegant approach to this?

L = 1.; (* step length *)
nMax = 1; (* number of realizations for averaging *)
nSteps = 1; (* number of modules in wire *)

(* define module *)
mod[x_] := UnitStep[x - L + 1.5]*UnitStep[L - x - 0.5]
(* define random profile *)
xi[x_] :=
Sum[RandomReal[{-0.5, 0.5}]*mod[x - n L], {n, -nSteps, nSteps}];
(* get nMax realizations of random profile *)
xiTmp[x_] = Table[xi[x], {n, 1, nMax}]
(* define correlator; normalized to corr[0]=1 *)
corr[y_]:=Table[Convolve[xiTmp[x][[n]], xiTmp[x][[n]], x, y]/
Convolve[xiTmp[x][[n]], xiTmp[x][[n]], x, 0], {n, 1, nMax}]
(* plot mean autocorrelator *)
Plot[Mean[corr[y]],{y,-5,5},PlotRange->All,Filling->Axis]

-
I'm not following your code but the slowdown seems be because of 0 in Convolve[xiTmp[x][[n]], xiTmp[x][[n]], x, 0] –  Mr.Wizard Jul 29 '12 at 18:04
Welcome to Mathematica SE @nilfisque! Why all these Table that run from 1 to 1? If there is a slowdown and you want to get to the root of it, simplifying the code to bare minimum for readability and easier debugging is appreciated. –  Vitaliy Kaurov Jul 29 '12 at 18:10
Could you elaborate? This Convolve[f,f,x,y]/Convolve[f,f,x,0] construct ensures that corr[y=0] is equal to unity and this is an essential requirement in my analysis. –  nilfisque Jul 29 '12 at 18:16
@VitaliyKaurov Thank you. I actually tried to scale down the code but apparently did a very bad job at it. The Table functions stem from my original calculations where I wanted to average over >100 realizations. These calculations didn't finish in a reasonable amount of time, so I set them to 1 for optimization. –  nilfisque Jul 29 '12 at 18:25
btw if you want to calculate autocorrelator of f[x], you may want to use Convolve[f[x], f[-x], x, y] instead of Convolve[f[x], f[x], x, y]. –  Silvia Jul 29 '12 at 19:36

First I want to say, as you mentioned in your comment that your ultimate goal is to to do it for nMax over 100, I suggest you first symbolicly calculate the correlation of the following function, treating $r_n$ ($n=-s,-s+1,\dots,s$, and $s$ is nSteps for short) as variables as $x$:

$$\xi(x,r_{-s},r_{-s+1},\dots,r_{s})=\sum _{n=-s}^{s} r_n\, \text{mod}(x-n\,L)\text{.}$$

On the othe hand, if what you need is only the final plot, maybe it's more convenient to do it in discrete style. Supplying each $\xi(x)$ as a discrete List, you can use numerical ListCorrelate on your data for great speed-up, avoiding symbolic calculation of Convolve for functions. Here we calculate the plot for nMax = 1000:

L = 1.; (* step length *)
nMax = 1000; (* number of realizations for averaging *)
nSteps = 1; (* number of modules in wire *)

(* define module *)
mod[x_] := UnitStep[x - L + 1.5]*UnitStep[L - x - 0.5]
(* define random profile *)
xi[x_] :=
Sum[RandomReal[{-0.5, 0.5}]*mod[x - n L], {n, -nSteps, nSteps}];
(* get nMax realizations of random profile *)
xiTmp[x_] = Table[xi[x], {n, 1, nMax}]


Now we discrete the xiTmp on a x-grid:

dx = .01;

{xiSupportInf, xiSupportSup} =
{(-(nSteps - 1))*L - 3/2, (nSteps + 1)*L - 1/2};

yWidth = xiSupportSup - xiSupportInf;

AbsoluteTiming[xiTmpDiscrete =
Transpose[Table[xiTmp[x], {x, xiSupportInf,
xiSupportSup, dx}]]; ]


{4.2972458, Null}

Lx = Dimensions[xiTmpDiscrete][[2]];

AbsoluteTiming[
corrArray = ((#1/#1[[Lx]] & )[ListCorrelate[#1,
#1, {-1, 1}, 0]] & ) /@ xiTmpDiscrete; ]


{0.0950055, Null}

ListLinePlot[
Transpose[
{Range[-yWidth, yWidth, dx],
Mean[corrArray]}
], PlotRange -> All,
Filling -> Axis]


This may still be slow for large nSteps, in that case you may want to somehow Compile the xiTemp function.

Update: some note about the correlate calculation

Here xiTmpDiscrete is a List whose elements are those $\xi(x)$s (in discrete form, so each $\xi(x)$ itself is a List too, in the form $\{\xi(x_1), \xi(x_2), \dots, \xi(x_{\mathrm{Lx}})\}$).

The pure function ListCorrelate[#1, #1, {-1, 1}, 0]& takes a $\xi(x)$ as its input, calculates the discrete version correlate $\int_{-\infty}^{\infty}\xi(x)\xi(x-y)\,\mathrm{d}x$.

Then the result is supplied to function #1/#1[[Lx]] &. As the ListCorrelate result with above setting is symmetric, the middle element, which is the Lx-th one, correspond to auto-correlate $\int_{-\infty}^{\infty}\xi(x)\xi(x-0)\,\mathrm{d}x$, so #1/#1[[Lx]] & effectively performs the normalization. Thus now we have obtained the normalized correlator for one $\xi(x)$.

By taking advantage of /@, which is a short-form for Map, the above process is performed for each $\xi(x)$ in xiTmpDiscrete. The final result corrArray is a array of size nMax $\times$ Lx. By applying Mean on it, mean is performed among rows, obtain a 1-dimension List of length Lx, which is a discrete version of your Mean[normalisedcorr[y]] on interval $[$ -yWidth, yWidth $]$.

-
Silvia, this is indeed way faster than my initial method, thanks a lot! May I ask you a further question about your construction of corrArray? While this line is very compact, I'm not sure I entirely understand it. So at first you apply the list xiTmpDiscrete to ListCorrelate and than the resulting expression is applied to #1/#1[[Lx]]. I've tried to cast this into code I'm more familiar with (at least right now) and came up with ListCorrelate[xiTmpDiscrete, xiTmpDiscrete, {-1, 1}, 0]/ ListCorrelate[xiTmpDiscrete, xiTmpDiscrete, {-1, 1}, 0][[Lx]], which however results in an error. –  nilfisque Jul 30 '12 at 15:12
@nilfisque please see my update. (it's too long for a comment.) –  Silvia Jul 30 '12 at 16:29

Here are two things you can do to speed up this code.

## 1. Do the convolution with symbolic y

Because you have defined corr using SetDelayed, the table of Convolve expressions will be re-evaluated every time you evaluate corr[number]. The normalisation term with y=0 is causing a particular slow down, though I'm not sure why exactly.

If you instead use Set the convolution will be computed straight away with symbolic y, which means that every time you evaluate corr[number] Mathematica will just take the expression it has already computed and insert number in place of y.

corr[y_] = Table[Convolve[xiTmp[x][[n]], xiTmp[x][[n]], x, y], {n, 1, nMax}];


You'll notice that I've removed the normalisation term. Now that we have a symbolic expression for corr[y] we can get the normalisation term simply by evaluating corr[0]. The normalised correlation can therefore be defined as:

normalisedcorr[y_] := corr[y]/corr[0]


## 2. Convert xi to a Piecewise expression

The convolution integrals can be made faster by expressing xi[x] as a Piecewise expression instead of in terms of UnitStep. To do this we can insert PiecewiseExpand into the definition of xi:

xi[x_] := PiecewiseExpand @ Sum[RandomReal[{-0.5, 0.5}]*mod[x - n L], {n, -nSteps, nSteps}]


For convenience, the whole tweaked code is:

L=1.;nMax=1;nSteps=1;
mod[x_]:=UnitStep[x-L+1.5]*UnitStep[L-x-0.5];
xi[x_]:=PiecewiseExpand@Sum[RandomReal[{-0.5,0.5}]*mod[x-n L],{n,-nSteps,nSteps}];
xiTmp[x_]=Table[xi[x],{n,1,nMax}];
corr[y_]=Table[Convolve[xiTmp[x][[n]],xiTmp[x][[n]],x,y],{n,1,nMax}];
normalisedcorr[y_]:=corr[y]/corr[0]
Plot[Mean[normalisedcorr[y]],{y,-5,5},PlotRange->All,Filling->Axis]


This takes about 19 seconds on my PC with nMax=100

-
Does this do the same thing? You completely removed the troublesome Convolve[xiTmp[x][[n]], xiTmp[x][[n]], x, 0]. Also, since the OP is using RandomReal I assumed caching was verboten. –  Mr.Wizard Jul 29 '12 at 18:58
@Mr.Wizard, I think so. The troublesome piece simply computes the convolution evaluated at y=0. Since we have already computed the convolution for all values of y, it's not really required. We can just normalise with corr[y]/corr[0] –  Simon Woods Jul 29 '12 at 19:05
@SimonWoods, thank you for your input on my code, I appreciate it a lot! These suggestions are exactly what I was hoping for (I'm especially excited about the elegant corr[y]/corr[0] construction. :)). –  nilfisque Jul 30 '12 at 15:34