Take the 2-minute tour ×
Mathematica Stack Exchange is a question and answer site for users of Mathematica. It's 100% free, no registration required.

I'm working on some pretty intense computation in Mathematica; when my code started running slowly, I tracked the source of the problem to Exp[]. I need to exponentiate every element of a 50x500x500 array; performing the operation on a 500x500 array takes on the order of 3 seconds (according to AbsoluteTime), so the entire array should take about 50 times that. Unfortunately, that's calculation needs to happen for every data point.

I've read about lots of ways to speed up Mathematica code, but none of those methods seem to apply here. I'm already working in MachinePrecision. I have noticed that some of my results are ridiculously small (for example, 4.282835067271648*10^-78127094), but I'm not sure how to make Mathematica ignore those; they're obviously much smaller than $MachineEpsilon.

Any advice is greatly appreciated!

Update:

Below is a sample of my code and the generated output. To give it some context, g0, is a scalar, σg0 is a length 50 array, and g is a 500x500 array.

(* Added after Oleksandr R.'s comment *)
SetSystemOptions["CatchMachineUnderflow" -> False];

n = Length[σg0];
probgs = ConstantArray[N[0], {50, 500, 500}];
For[i = 1, i <= n, i++,
  probgs[[i]] = 
    N[(1/(Sqrt[2 π] σg0[[i]])) Exp[-0.5 ((g - g0)/σg0[[i]])^2]];
  ]; // AbsoluteTiming
Precision[probgs]

Output:

{4.816275, Null}
MachinePrecision

Turning off underflow definitely helped; 5 seconds isn't bad at all for what I'm doing.

share|improve this question
1  
xxx = RandomReal[{-10, 10}, {50, 500, 500}]; Exp[xxx]; takes about 0.04 seconds for me. –  b.gatessucks May 8 '13 at 10:42
7  
4.282835067271648*10^-78127094 is not a MachinePrecision number, so actually you are working in arbitrary precision. You have been caught out by automatic underflow handling, whereas clearly you would rather that such numbers go to (machine) zero. Evaluate SetSystemOptions["CatchMachineUnderflow" -> False] and try again. –  Oleksandr R. May 8 '13 at 10:50
4  
If you're evaluating Exp[] at very tiny arguments, you might consider replacing Exp[] with a suitably chosen Padé approximant. –  J. M. May 8 '13 at 11:22
3  
When dealing with probabilities, sometimes it can be more efficient to work with their natural logarithm instead of the probabilities themselves. –  Thies Heidecke May 8 '13 at 11:33
2  
What is \[Sigma]g0? Standalone means it should run from a fresh kernel w/o any alterations. The actual arrays matter and you cannot expect other users to bother to generate dummys. –  Yves Klett May 8 '13 at 12:45

3 Answers 3

up vote 14 down vote accepted

Obviously, for large negative inputs, Exp will produce very small numbers. While this isn't intrinsically problematic, it so happens that, by default, Mathematica deals with machine underflow by converting the affected values to an arbitrary precision representation in order to avoid catastrophic loss of precision. However, sometimes one would rather disregard underflowed values instead (i.e. let them go to zero), and indeed that seems to be the case here.

This behavior can be controlled using the system option "CatchMachineUnderflow"--simply use

SetSystemOptions["CatchMachineUnderflow" -> False]

and underflowed values will be flushed to (machine precision) zero.

Since this is a global option that will most likely affect the results of system functions as well as user code, it's advisable to localize its effect as tightly as possible. For this purpose one can use the undocumented function Internal`WithLocalSettings, as described by Daniel Lichtblau in this StackOverflow answer:

With[{cmuopt = SystemOptions["CatchMachineUnderflow"]},
 Internal`WithLocalSettings[
  SetSystemOptions["CatchMachineUnderflow" -> False],
  (* put your own code here; for example: *)
  Exp[-1000.],
  SetSystemOptions[cmuopt]
 ]
]
(* 0.` *)

Contrast this with:

Exp[-1000.]
(* 5.0759588975494567652918094795743369258164499728`12.954589770191006*^-435 *)
share|improve this answer

Unfortunately, my Tarot-Cards are in repair so here is a guess:

data = yourSecertData;
Precisionp[data]
pd = Developer`ToPackedArray[N[data, $MachinePrecision]];
(* must be True *)
Developer`PackedArrayQ[pd]
AbsoluteTiming[pd];
cExp = Compile[{{in, _Real, 2}}, Exp[in], 
   RuntimeAttributes -> Listable, CompilationTarget -> "C"];
AbsoluteTiming[cExp[pd]]

You may also want to play with the "Listable" setting and the tensor rank (2) and see if parallel execution helps, or if compilation to "C" is needed. (Maybe "WVM" is enough). If using the compiler seems helpful, may want to look up other compiler options (RuntimeOptions).

Update:

After more info was available:

cf = With[{g0 = g0}, 
  Compile[{{og0, _Real, 0}, {g, _Real, 
     2}}, (1/(Sqrt[2 \[Pi]] og0)) Exp[-0.5 ((g - g0)/og0)^2], 
   RuntimeAttributes -> Listable]]

This should speed up things considerably.

share|improve this answer
5  
Compiling can be just as esoteric as the tarot cards, though :-) –  Yves Klett May 8 '13 at 10:51
    
@YvesKlett, :-) lol. –  user21 May 8 '13 at 10:53
    
@ruebenko, I tried messing around with Compile for a while, but it complained about a pure function I had defined and, in the end, crashed my system. I never figured out what I was doing wrong, unfortunately. –  nosuchthingasstars May 8 '13 at 11:14
    
@nosuchthingasstars, it may actually be helpful if you post that code such that this can be verified/fixed/improved. None can answer you question without seeing the code. This is like saying my bike is slow, can fix it? –  user21 May 8 '13 at 11:17
    
@ruebenko, I added the code just now. What I posted should run standalone. –  nosuchthingasstars May 8 '13 at 11:20

You will save some time by computing -0.5 (g-g0)^2 just once:

probgs = With[{gg = -0.5 (g - g0)^2},
    Table[1/(Sqrt[2 Pi] sg0[[i]]) Exp[gg / sg0[[i]]^2], {i, n}]];
share|improve this answer

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.

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