5
$\begingroup$

In contrast to earlier versions, Mathematica 11.3 and 12.0 generate lots of underflow errors from my code, e.g. 

"General::munfl: Exp[-11250.] is too small to represent as a normalized machine number; precision may be lost."

A minimal working example: 

sigma=0.01;
gaussian=Function[x,1/(Sqrt[2Pi] sigma) Exp[-1/2(x/sigma)^2]]
gaussianP[x_]=D[gaussian[x],x]
gaussianP[1.5]

The background of this error is discussed in:

New General::munfl error and loss of precision How to flush machine underflows to zero and prevent conversion to arbitrary precision?

In short: Mathematica still detects MachineNumber underflows, but it does no longer switch to arbitrary precision mode. Instead, it produces a result with less digits if any.

The solution proposed in question 170416 does not solve my problem since Exp[x] expressions are converted to E^x upon evaluation, cf. https://reference.wolfram.com/language/ref/Exp.html. Thus it is not sufficient to only modify Exp[].

Question 1: How can I restore the behavior of previous versions?

Considering the solution in 170416, how can I tell Mathematica to treat E^MachineNumber in the same way as Exp[MachineNumber]?

There are already attempts to rewrite E^x expressions in terms of Exp[x]: Format Exp[] in output Function to Expand exponentials

How can these be used in my case?

Question 2: How can I handle underflow manually?

I feel uncomfortable with Mathematica producing results of unknown precision. Instead, I would like to define locally a threshold by myself, below which the outcome of Exp[] is set to zero at MachinePrecision.

I would be grateful for any help.

Improve this question
$\endgroup$
6
  • $\begingroup$ In general any error or warning can handled manually with Check $\endgroup$
    – m_goldberg
    May 5, 2019 at 20:14
  • 1
    $\begingroup$ related: mathematica.stackexchange.com/q/197580/26598 $\endgroup$
    – Roman
    May 5, 2019 at 20:28
  • 1
    $\begingroup$ I think that Mathematica now underflows in accordance with IEEE 754. Can anyone confirm? $\endgroup$
    – Michael E2
    May 5, 2019 at 22:13
  • $\begingroup$ The point of my comment is that Mathematica does not produce "results of unknown precision," but one that conforms to a defined standard. However, I have not found any documentation that explicitly confirms. This tutorial describes the new underflow behavior in a way that seems consistent with IEEE 754. $\endgroup$
    – Michael E2
    May 6, 2019 at 11:55
  • 1
    $\begingroup$ @Michael E2 much of that tutorial is agnostic about the float point system used, and describes the constants that may be employed to characterize the system. But the part that describes the new underflow behavior seems specific to IEEE 754, since denormalized numbers are a special feature of that particular standard. Perhaps Mathematica no longer supports other floating point systems. $\endgroup$
    – John Doty
    May 6, 2019 at 13:00

2 Answers 2

Reset to default
4
$\begingroup$

The old behavior was inconsistent with the principle that machine numbers have the semantics of your hardware's floating point. The new behavior conforms to this principle. Thus, if you want controlled precision in your calculations, you must use controlled precision as input, or explicitly set the precision at strategic places.

Improve this answer
$\endgroup$
5
  • 1
    $\begingroup$ BTW, it's still inconsistent with respect to overflow, though. Compare 10.^-400 and 1/10.^400. $\endgroup$
    – Michael E2
    May 6, 2019 at 11:49
  • 2
    $\begingroup$ @MichaelE2 That's because there is no "obvious" machine number to overflow to. NaN might be possible, but overflowing to a bignum tends to be more useful (maybe you knew this, but it perhaps is not obvious to others). $\endgroup$
    – Daniel Lichtblau
    May 6, 2019 at 16:58
  • 1
    $\begingroup$ @DanielLichtblau IEEE specifies overflow should round to ± infinity or the maximum machine number, depending on rounding mode. Round-to-even would have it be (floating-point) infinity. But (I suppose) Mathematica's Infinity is not equivalent to machine floating-point infinity. Something like an equivalent MachineInfinity or DirectedInfinity[1.] would have to be introduced and incorporated into WL handle machine overflow, I guess. (It'd be nice if packed arrays could contain such floating-point numbers.) $\endgroup$
    – Michael E2
    May 6, 2019 at 18:30
  • $\begingroup$ Also I agree that overflowing to bignums is useful. I thought that about underflowing, too. $\endgroup$
    – Michael E2
    May 6, 2019 at 18:31
  • 2
    $\begingroup$ @MichaelE2 Right, DirectedInfinity[1] is not equivalent to machine float infinity. Possibly a float infinity for Mathematica would be useful though. As for underflowing to a bignum, the pains outweigh the gains. $\endgroup$
    – Daniel Lichtblau
    May 6, 2019 at 20:47
4
$\begingroup$

One idea is to use Check to find out when General::munfl happens, and then switch to a more complicated evaluation process. The more complicated evaluation process would modify the evaluations of numeric functions so that they rationalize their input first, and then numericize the output. Here is a function that does this:

SetAttributes[trappedEvaluation, HoldFirst]

trappedEvaluation[expr_] := Check[
    expr,
    numericize @ Block[{Real},
        f_?nfQ[a___, r_Real, b___] ^:= Rationalize[Unevaluated @ f[a, r, b], 0];
        expr
    ],
    General::munfl
]

nfQ = System`Dump`HeldNumericFunctionQ;

numericize[e_] := With[{res = N[e, 16]},
    Quiet @ Check[
        N[res],
        res,
        General::munfl
    ]
]

Your example:

trappedEvaluation[gaussianP[1.5]]

General::munfl: Exp[-11250.] is too small to represent as a normalized machine number; precision may be lost.

-9.206189566469410*10^-4881

Improve this answer
$\endgroup$
5
  • $\begingroup$ I think you forgot a SetAttributes[trappedEvaluation, HoldAll] or something like it. $\endgroup$
    – Michael E2
    May 6, 2019 at 20:03
  • $\begingroup$ @MichaelE2 Good point! thanks $\endgroup$
    – Carl Woll
    May 6, 2019 at 20:05
  • $\begingroup$ Very cool, this helped me a lot! Is it also possible to use the rationalize-numericize trick when storing an intermediate result in a symbol, something like val = trappedEvaluation[gaussianP[1.5]]; trappedEvaluation[2. val]? $\endgroup$
    – Quantum Sharpener
    May 8, 2019 at 14:37
  • $\begingroup$ Let's assume I want to calculate res = trappedEvaluation[2. gaussianP[1.5]] in two steps: 1) val = trappedEvaluation[gaussianP[1.5]] (*val is of arbitrary precision*) 2) res = 2. val This results in 0. (machine precision!) for res associated with an underflow error. How can I prevent that? $\endgroup$
    – Quantum Sharpener
    May 14, 2019 at 21:44
  • $\begingroup$ @CarlWoll Thank you; I've found this very helpful. However, when I evaluate trappedEvaluation[N[Exp[-800]]], it still returns zero. However, if I define f[x_?NumberQ] = N[Exp[-x]], and then run trappedEvaluation[f[800]] it returns the correct value. Could you please explain what goes wrong in the first case? $\endgroup$
    – Lauren Pearce
    Mar 12, 2021 at 0:09

Your Answer

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

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