# Efficient evaluation of function defined by multiple Sum

I tried to make a workable example here. The problem is that it is very very slow. I must be doing something very stupid out here. Can somebody give me some clue how to evaluate this multiple sum faster. Right now it takes around 2 mins. However, for real problem I need the parameter Npol to be of the order of ~10.

Clear["Global*"];
SetSystemOptions["CompileOptions"->"CompileReportExternal"->False];
DRSH=N[(94688879/10000000)*10^(-7)];
Npol=2;
nk=40;
PN=Evaluate[Table[0.0025 n,{n,1,Npol}]];
OmN=Evaluate[Table[(1.0/DRSH)^((-n+2)/Npol),{n,1,Npol}]];
fpol[n_, On_]:=(-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n)/; On== OmN[[1]]
fpol[n_, On_]:=(-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n) (0.77 2.0^-n-On)/(0.77 2.0^-n+On)

OverLaps1=Compile[{{k, _Integer},{n,_Real}, {m,_Real}},Exp[-0.5 Sum[If[q!= k,(fpol[q, n]-fpol[q, m])^2,0],{q,0,nk}]],CompilationTarget->"C"];
OverLaps2=Compile[{{k, _Integer},{n,_Real}, {m,_Real}},Exp[-0.5 Sum[If[q!= k,(fpol[q, n]+fpol[q, m])^2,0],{q,0,nk}]],CompilationTarget->"C"];
NpolSum=Compile[{{i, _Integer},{j, _Integer}, {k, _Integer}}, (Sum[ Sum[PN[[n]]  PN[[m]] fpol[k, OmN[[n]]]^i fpol[k, OmN[[m]]]^j Exp[-0.5 (fpol[k, OmN[[n]]]^2+ fpol[k, OmN[[m]]]^2)] ( OverLaps1[k, OmN[[n]], OmN[[m]]]+ (-1.0)^j   OverLaps2[k, OmN[[n]], OmN[[m]]]),{n,1,Npol}],{m,1,Npol}])^2,CompilationTarget->"C"];

S2Entrp=Compile[{{k, _Integer}},Sum[Sum[ ( 1/(i! j!)  )  NpolSum[i, j, k],{i,0,10}],{j,0,10}]];
ListLogLinearPlot[Table[{k,S2Entrp[ k]}, {k,0, nk}],PlotStyle->Directive[PointSize[Medium],Red], PlotRange-> All]//AbsoluteTiming


I thank you for your time.

-
Without looking at the code in detail, did you try replacing Sum[...,{}] by Total@Table[...,{}]? – acl Mar 29 '13 at 18:59

One significant thing you could do is memoize the definitions for fpol and both OverLaps functions, as these are called multiple times with the same arguments.

I would also suggest creating a definition directly for fpol[n_, OmN[[1]]] rather than using the condition on the RHS.

I have ditched the Compile here, as I don't think it's giving any benefit due to all the calls to MainEvaluate (BTW suppressing the warnings doesn't make the problem go away...)

fpol[n_, OmN[[1]]] := fpol[n, OmN[[1]]] = (-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n)
fpol[n_, On_] := fpol[n, On] = (-0.38 Sqrt[4^-n])/(DRSH + 0.77 2.0^-n) (0.77 2.0^-n - On)/(0.77 2.0^-n + On)

OverLaps1[k_, n_, m_] := OverLaps1[k, n, m] = Exp[-0.5 Sum[If[q != k, (fpol[q, n] - fpol[q, m])^2, 0], {q, 0, nk}]]
OverLaps2[k_, n_, m_] := OverLaps2[k, n, m] = Exp[-0.5 Sum[If[q != k, (fpol[q, n] + fpol[q, m])^2, 0], {q, 0, nk}]]

NpolSum[i_, j_, k_] := (Sum[Sum[PN[[n]] PN[[m]] fpol[k, OmN[[n]]]^i fpol[k, OmN[[m]]]^j
Exp[-0.5 (fpol[k, OmN[[n]]]^2 + fpol[k, OmN[[m]]]^2)] (OverLaps1[k, OmN[[n]],
OmN[[m]]] + (-1.0)^j OverLaps2[k, OmN[[n]], OmN[[m]]]), {n, 1, Npol}], {m, 1, Npol}])^2

S2Entrp[k_] := Sum[Sum[(1/(i! j!)) NpolSum[i, j, k], {i, 0, 10}], {j, 0, 10}]


I think there are further improvements that could be made, for example making fpol and OverLaps Listable and replacing NpolSum with something nifty using Outer, but the memoization appears to be the key to increasing the performance here.

With these changes the plot takes under half a second to create:

ListLogLinearPlot[Table[{k, S2Entrp[k]}, {k, 0, nk}],
PlotStyle -> Directive[PointSize[Medium], Red], PlotRange -> All] // AbsoluteTiming
`

-
Yes,thank you. I also found out that Compile doesn't help here. Again thank you for the quick tip. – SSB Mar 29 '13 at 23:43