1
$\begingroup$

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.

$\endgroup$
1
  • $\begingroup$ Without looking at the code in detail, did you try replacing Sum[...,{}] by Total@Table[...,{}]? $\endgroup$
    – acl
    Mar 29, 2013 at 18:59

1 Answer 1

1
$\begingroup$

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

enter image description here

$\endgroup$
1
  • $\begingroup$ Yes,thank you. I also found out that Compile doesn't help here. Again thank you for the quick tip. $\endgroup$
    – SSB
    Mar 29, 2013 at 23:43

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.