Skip to main content

performance tuning - How to flush machine underflows to zero and prevent conversion to arbitrary precision?


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.



Answer




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 *)

Comments

Popular posts from this blog

plotting - How to draw lines between specified dots on ListPlot?

I would like to create a plot where I have unconnected dots and some connected. So far, I have figured out how to draw the dots. My code is the following: ListPlot[{{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4,13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full] I have thought using ListLinePlot command, but I don't know how to specify to the command to draw only selected lines between the dots. Do have any suggestions/hints on how to do that? Thank you. Answer One possibility would be to use Epilog with Line : ListPlot[ {{1, 1}, {2, 2}, {3, 3}, {4, 4}, {1, 4}, {2, 5}, {3, 6}, {4, 7}, {1, 7}, {2, 8}, {3, 9}, {4, 10}, {1, 10}, {2, 11}, {3, 12}, {4, 13}, {2.5, 7}}, Ticks -> {{1, 2, 3, 4}, None}, AxesStyle -> Thin, TicksStyle -> Directive[Black, Bold, 12], Mesh -> Full, Epilog -> { Line[ ...

dynamic - How can I make a clickable ArrayPlot that returns input?

I would like to create a dynamic ArrayPlot so that the rectangles, when clicked, provide the input. Can I use ArrayPlot for this? Or is there something else I should have to use? Answer ArrayPlot is much more than just a simple array like Grid : it represents a ranged 2D dataset, and its visualization can be finetuned by options like DataReversed and DataRange . These features make it quite complicated to reproduce the same layout and order with Grid . Here I offer AnnotatedArrayPlot which comes in handy when your dataset is more than just a flat 2D array. The dynamic interface allows highlighting individual cells and possibly interacting with them. AnnotatedArrayPlot works the same way as ArrayPlot and accepts the same options plus Enabled , HighlightCoordinates , HighlightStyle and HighlightElementFunction . data = {{Missing["HasSomeMoreData"], GrayLevel[ 1], {RGBColor[0, 1, 1], RGBColor[0, 0, 1], GrayLevel[1]}, RGBColor[0, 1, 0]}, {GrayLevel[0], GrayLevel...

list manipulation - Selecting multiple columns from a matrix?

Sample data: data = { {{2013, 1, 1}, 24.13, 167.67, 231.82}, {{2013, 1, 2}, 32.15, 170.92, 225.99}, {{2013, 1, 3}, 35.43, 172.68, 221.67}, {{2013, 1, 4}, 36.73, 173.05, 218.32}, {{2013, 1, 5}, 58.19, 165.96, 197.05}, {{2013, 1, 6}, 69.99, 163.50, 187.52}, {{2013, 1, 7}, 71.37, 154.21, 175.58}, {{2013, 1, 8}, 72.51, 149.66, 163.25}}; I want a DateListPlot with three graphs, so for a matrix formed by columns 1 and 2, one for columns 1 and 3, and 1 for columns 1 and 4. At the moment I'm using this code: data2 = Transpose[{data[[All, 1]], data[[All, 2]]}]; data3 = Transpose[{data[[All, 1]], data[[All, 3]]}]; data4 = Transpose[{data[[All, 1]], data[[All, 4]]}]; DateListPlot[{data2, data3, data4}, Joined -> True, Filling -> {3 -> {1}}] but I have a hunch that this can be done more efficiently. I don't like the Transpose s in particular. Any ideas? edit (for extra credit) What if I need to multiply the second column by 2, which in my solution is simp...