5
$\begingroup$

This question,which is still unanswered, might be relevant because it involves NIntegrate over lists and it also has Exp.

In the question linked I encountered slow evaluation of NIntegrate. I thought it might be because Mathematica somehow does the calculations slowly over lists and it is not the fault of NIntegrate the slow performance, however I am not sure yet. Anyway, I tried to calculate the integral by the simplest method that is boxing the area under integral. Consider this:

<<Developer`
axis = ToPackedArray@Table[1.*i, {i, 0, 10000, 0.01}];
Total[0.01*Exp[-axis]] // AbsoluteTiming

{10.108011, 1.0050083333194446}

It is really slow compared to Matlab. Then I tried:

Exp[-axis]; // AbsoluteTiming

{8.998900, Null}

and

Sin[-axis]; // AbsoluteTiming

{0.009001, Null}

It seems exponential on a list is really slow.

Compiling helps a lot:

exp = Compile[{{axis, _Real, 1}},

  Apply[Plus, 0.01*Exp[-axis]],
  Parallelization -> True, CompilationTarget -> "C"
  ]

exp[axis] // AbsoluteTiming

{0.062006, 1.0050083333194315}

But still Matlab is ten time faster.

Matlabs code:

axis = 0:0.01:10000;
tic
sum(exp(-axis).*0.01)
toc

ans =

1.0050

Elapsed time is 0.006555 seconds.

Why is Mathematica slow in this case and how can I make it faster?

I think this happens because Mathematica keeps really small numbers:

Exp[-1000000.]

3.296831478*10^-434295

How can I prevent this? Is it related to the issue at all?

Edit

It seems for Mathematica length isn't important but size does matter!

axis2 = Table[Abs@Sin[1.*i], {i, 0, 10000, 0.01}];

Total[0.01*Exp[-axis2]] // AbsoluteTiming

{0.015000, 5558.330252144095}

Although the length of the list is the same as before its values are between zero and one and it's faster.

$\endgroup$
  • $\begingroup$ Interesting question. I don't think Packeing will help in this case because Exp breaks this Packing. PackedArrayQ[Exp[-axis]] is False. $\endgroup$ – Algohi Mar 17 '15 at 15:34
  • $\begingroup$ Possibly related: Exp is slower for some number ranges $\endgroup$ – Karsten 7. Mar 17 '15 at 18:17
  • 1
    $\begingroup$ Matlab returns 0 for exp(-1000000.). $\endgroup$ – Karsten 7. Mar 17 '15 at 18:25
  • 3
    $\begingroup$ Related: "CatchMachineUnderflow" -> False. $\endgroup$ – Michael E2 Mar 17 '15 at 21:16
  • $\begingroup$ @Karsten7 I didn't get the point. Did you mean the implementation of Exp in Mathematica is less efficient? $\endgroup$ – MOON Mar 17 '15 at 21:53
6
$\begingroup$

Starting with an extension to the point raised in your edit.

rL = RandomReal[{0.5, 1.5}, 10000];
timings = Table[{10^n, AbsoluteTiming[Exp[-(rL*10^n)]][[1]]}, {n, 0, 5, .01}];
ListLogLinearPlot[{timings, {{746, 0}, {746, 0.05}}}, 
  AxesLabel -> {"Scaling", "Time"}, Joined -> {False, True}, PlotStyle -> {Black, Red}]

Plot

The red vertical line is placed at n=746. Matlab returns 0. for exp(-n) if n>=746.

Setting

SetSystemOptions["CatchMachineUnderflow" -> False]

as shown in this answer by Oleksandr R., for which Michael E2 provided a link in his comment to this question, improves the timings significantly.

Plot2

Mathematica is slower than Matlab, because Matlab uses machine precision whereas in Mathematica Exp can be evaluated to arbitrary numerical precision.

One can define a numerical version of Exp to "make it faster" using

nExp = Compile[{{x, _Real, 0}},
        Exp[x], CompilationTarget -> "C", Parallelization -> True, 
        RuntimeAttributes -> Listable, RuntimeOptions -> "Speed"]

Now

{nExp[-745], nExp[-746]}

{5.*10^-324, 0.}

similar to Matlab, the timings are independent of the scaling numericExp

and

0.01 Total[nExp[-axis]] // Timing

{0., 1.00501}

is much faster.

$\endgroup$
  • $\begingroup$ BTW, NSum[0.01 Exp[-i], {i, 0., 10000., 0.01}] // Timing and Sum[0.01 Exp[-i], {i, 0., 10000., 0.01}, Method -> "PolynomialExponential"] // Timing are both returning {0., 1.00501}. $\endgroup$ – Karsten 7. Mar 18 '15 at 1:22

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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