Systematically Optimizing Code for Better Runtime

Can someone please tell me how can I learn to systematically optimize my Mathematica code to be faster? I have seen a blog about the same but have no idea where to begin and how to implement this.

For example :- The following code takes more than 24 hours to run (actually it's still running, expected 30+ hours)

Ef[a_] := Pi^2*(a + 2)^2;
Eb[a_] := Pi^2*(a + 1)^2;
f[n_, x_] := Sqrt[2/((n + 2)^2 - 1)]*((n + 2)*Cos[Pi*(n + 2)*x] - Cot[Pi*x]*Sin[(n + 2)*Pi*x]);
b[n_, x_] := Sqrt[2]*Sin[(n + 1)*Pi*x];
xf[m_, n_] := If[Mod[m - n, 2] == 0, 0, Integrate[f[n, x]*f[m, x]*x, {x, 0, 1},
Assumptions -> {Element[n, Integers], Element[m, Integers]}]];
xb[m_, n_] := Integrate[b[n, x]*b[m, x]*x, {x, 0, 1}, Assumptions -> {Element[n, Integers], Element[m, Integers]}];
xt[m_, n_] := If[m == 0 || n == 0, xb[m, n], (1/2)*(xb[m, n] + xf[m - 1, n - 1])];
Z[T_] := Sum[E^(-(Eb[i]/T)), {i, 0, 10}];
y1[m_, t_] := -Sum[(Eb[k] - Eb[m])*xt[m, k]*xt[k, m]*Cos[(Eb[k] - Eb[m])*t], {k, 0, 10}];
Y1[T_, t_] := (-(1/Z[T]))*Sum[Sum[(Eb[k] - Eb[m])*xt[m, k]*xt[k, m]*Cos[(Eb[k] - Eb[m])*t], {k, 0, 10}]/E^(Eb[m]/T),
{m, 0, 10}];
file = OpenAppend["susypotwell_Y1T0.1.dat"];
Table[Export[file, {{t, Y1[0.1, t]}}, "TSV"], {t, -1., 1., 0.01}]
Close[file]
file = OpenAppend["susypotwell_Y1T1.dat"];
Table[Export[file, {{t, Y1[1, t]}}, "TSV"], {t, -1., 1., 0.01}]
Close[file]
file = OpenAppend["susypotwell_Y1T10.dat"];
Table[Export[file, {{t, Y1[10, t]}}, "TSV"], {t, -1., 1., 0.01}]
Close[file]


One way to do this could be to evaluate Y1[T,t] and define a new function using the output which would then be used in the Export expressions instead of Y1. Please note that here I am exporting data to plot because Mathematica does not automatically save data of plots if it needs to be modified later on. If Plot can somehow do this faster then I have no problem with it as well and I will just use this or similar to save the plot data inside notebook for later manipulations.

Furthermore, if say, one has optimized their code as much as it can be and it is still taking 24+ hours or something like that then what are the options they have?

Apologies for such a broad question. Even quick tips and suggestions would help a lot.

Edit :-

Step-1 : Use #-& notation for pure functions. (source)

• “Exportin data to plot” - then I would strongly recommend turning to numerical calculations once you have substituted appropriate numerical values into your equations. After all, you wouldn’t be able to plot them without doing that anyway. Also, your idea of saving the result of Integrate and Sum once instead of recalculating it every time might be the most significant change you could make. – MarcoB Jun 24 at 1:58

Instead of 30 hours, do the job in 2.6 seconds.

Edit Used NIntegrate, as @flinty recommended.

    ClearAll["Global*"]
(Ef[a_] = Pi^2*(a + 2)^2;
Eb[a_] = Pi^2*(a + 1)^2;
f[n_, x_] =
Sqrt[2/((n + 2)^2 - 1)]*((n + 2)*Cos[Pi*(n + 2)*x] -
Cot[Pi*x]*Sin[(n + 2)*Pi*x]);
b[n_, x_] = Sqrt[2]*Sin[(n + 1)*Pi*x];

(*    Table[xf[m, n] =
If[Mod[m - n, 2] == 0, 0,
Integrate[f[n, x]*f[m, x]*x, {x, 0, 1}]], {m, 0, 10}, {n, 0, 10}];
Table[xb[m, n] = Integrate[b[n, x]*b[m, x]*x, {x, 0, 1}], {m, 0,
10}, {n, 0, 10}];   *)

Table[xf[m, n] =
If[Mod[m - n, 2] == 0, 0,
NIntegrate[f[n, x]*f[m, x]*x, {x, 0, 1}, MaxRecursion -> 50]], {m,
0, 10}, {n, 0, 10}];
Table[xb[m, n] =
If[(1/2 (2 + m + n)) \[Element] Integers && m != n, 0,
NIntegrate[b[n, x]*b[m, x]*x, {x, 0, 1}]], {m, 0, 10}, {n, 0,
10}];

xt[m_, n_] =
If[m == 0 || n == 0, xb[m, n], (1/2)*(xb[m, n] + xf[m - 1, n - 1])];
Z[T_] = Sum[E^(-(Eb[i]/T)), {i, 0, 10}] // Simplify;
y1[m_, t_] = -Sum[(Eb[k] - Eb[m])*xt[m, k]*xt[k, m]*
Cos[(Eb[k] - Eb[m])*t], {k, 0, 10}];
Y1[T_, t_] := (-(1/Z[T]))*
Sum[Sum[(Eb[k] - Eb[m])*xt[m, k]*xt[k, m]*
Cos[(Eb[k] - Eb[m])*t], {k, 0, 10}]/E^(Eb[m]/T), {m, 0, 10}];
tab1 = Table[{t, Y1[0.1, t]}, {t, -1., 1., 0.01}];
tab2 = Table[{t, Y1[1, t]}, {t, -1., 1., 0.01}];

tab3 = Table[{t, Y1[10, t]}, {t, -1., 1., 0.01}];
{ListLinePlot[tab1, Epilog -> {Red, Point@tab1}],
ListLinePlot[tab2, Epilog -> {Red, Point@tab2}],
ListLinePlot[tab3, Epilog -> {Red, Point@tab3}]}
) // Timing


Let me say it repeadetly: Avoid SetDelayed (:=) whereever you can. My opinion.

• One could go further and replace the Table's with Scan to avoid storing the intermediate results, and replace Integrate's with NIntegrate's - but this is more work and 92s is pretty good so far. – flinty Jun 24 at 12:27
• @flinty, besides regarding the symmetry in Table, i don't see at the moment, how to take advantage of Scan. Could you please help me with example code. – Akku14 Jun 24 at 14:04
• Glad to see NIntegrate made it so much faster. There's an example in the documentation like this: Scan[(u[#] = x) &, {55, 11, 77, 88}] but I've noticed a problem with Scan when iterating over multidimensional data. I will ask about this in a separate question. – flinty Jun 24 at 14:18
• First let me say, the immideate calculation with Table und Sum, where all parts get stored, is only practicable for not too high m and n. You can also define Y1[T_,t_] = ... with Set. SetDelayed here was a relict of my testing. For NIntegrate, i set all integrals, that are zero, explicitly to zero, because NIntegrate complains of failed convergence with machine precision and needs much longer for that integrals. If you do xt[m, n]  you see the if condition is simply stored to be calculated later for definite m and n. – Akku14 Jun 24 at 19:02
• Table of xf and xb indeed stores 11*11 values. Notice xt[m,n] = ... is not a pattern definiton, but definition for the discrete m and n given by Table (or Sum). Mathematica finds it, because it is stored under for example xf[2,5]. The main bootleneck was the calculation of the integral with Integrate in your code since it was called very often and recalculated again. My approach with Table and Set indeed consumes more memory space, but looking up already calculated valures is much faster. – Akku14 Jun 24 at 19:12

In order to determine where to concentrate your efforts, you need to know where your bottlenecks are.

To do this work through a single Y1 calculation step by step. I would target your Integrate & Sum.

I'd be writing all results at once instead of using OpenAppend too.

As a guide, on my Linux 18.04 XUbuntu 12.0 combination (Xeon E5-2690 v4 @ 2.60GHz), the calculation Y1[10., 1.] takes 221.04 seconds.

For Y1[1., 1.] I get lots of underflows. As an example

Timing[Y1[1., #]]& /@ {-1, 0, 1}
During evaluation of In[20]:= General::munfl: Exp[-799.438] is too small to represent as a normalized machine number; precision may be lost.
During evaluation of In[20]:= General::munfl: Exp[-986.96] is too small to represent as a normalized machine number; precision may be lost.
During evaluation of In[20]:= General::munfl: Exp[-1194.22] is too small to represent as a normalized machine number; precision may be lost.
During evaluation of In[20]:= General::stop: Further output of General::munfl will be suppressed during this calculation.
Out[20]= {{226.588,0.247026},{224.858,-0.999189},{224.499,0.247026}}


so the Y1[1., #]& /@ Range[-1., 1, .01] calculation is going to take (at an average of 226 seconds per t) about 12 hours and 40 minutes.

Failing this, I reach for gfortran.

• OpenAppend` is actually important because if there is a power cut etc and the computation stops then the computed data is not lost. Also, can you say what you mean by gfortran usage? I have really long and complicated expressions (besides the one given here) and inputting them in gfortran etc would pose more of a problem than this one. – Nitin Jun 24 at 2:59
• @Nitin: by gfortran, I mean Fortran. Your preference for compiled language may differ. My implication is that despite its considerable power, MMA is not suitable for every problem, and yours may be one of those problems. – dwa Jun 24 at 3:29
• The expressions I need to use will take more time to put into Fortran and verify that there is no error compared to the time gained by preferring it over Mathematica. Any other option you can think of? – Nitin Jun 24 at 3:46
• @Nitin: as I suggest, work through your calculation step by step. You use Sum & Integrate. Can you get away with NSum and NIntegrate? Do you need each of the 201 steps? Can you get away with using fewer steps, possibly in conjunction with Interpolation? – dwa Jun 24 at 3:51
• If you are looking for a solution that is going to take substantially less than 30 hours to implement and run then I suspect there are few or no options. You might buy a computer that is more than 3x faster than what you have. Or buy a computer with lots of cpu cores and 16x the memory than you have now and learn all the tricky details of doing MMA parallel computations. But neither of those are going to be implemented and give you your solution in less than 30 hours total AND you will have to deal with the issue of moving MMA from one licensed computer to another. – Bill Jun 24 at 3:55