70
$\begingroup$

How to Control the Precision and Accuracy of Numerical Results
Arbitrary-Precision Numbers

Mathematica works with exact numbers and with two different types of approximate numbers: machine-precision numbers that take advantage of specialized hardware for fast arithmetic on your computer, and arbitrary-precision numbers that are correct to a specified number of digits.


To be sure of n correct digits, use N[expr, n].


When you do a computation, Mathematica keeps track of which digits in your result could be affected by unknown digits in your input. It sets the precision of your result so that no affected digits are ever included. This procedure ensures that all digits returned by Mathematica are correct, whatever the values of the unknown digits may be.


Mathematica automatically increases the precision that it uses internally in order to get the correct answer

Of course, this sounds very reassuring, but I still have some doubts that all decimal digits ever returned by Mathematica when working with arbitrary-precision numbers are always provably correct, no matter what functions I invoked.

What are those cases when I can be certainly sure all displayed digits are correct?


Update:
Here is an example when some incorrect decimal digits are returned when working with arbitrary-precision arithmetic:

a = 1`7
(* 1.000000 *)

a // Precision
(* 7. *)

d = Derivative[0, 1][StieltjesGamma][0, a]
(* -1.6450 *)

MachineNumberQ[d]
(* False *)

d // FullForm
(* 1.64501552391043694947251282378009083269`5.155856939311388 *)

d // Precision
(* 5.15586 *)

So, Mathematica claims that at least 5 (hence, all) decimal digits of the result -1.6450 are correct. But in fact, the exact result is -Pi^2 / 6 that is -1.644934..., so only 3 digit are correct.

I am also concerned that Precision[...] itself returns a machine-precision number, that is subject to uncontrolled error-accumulation that possibly can result in claiming more digits of precision in a number than there actually are. Can I assume that Mathematica always errs on the safe side when computing a precision?


Update 2:
Another (gross) example:

a = 2`6
(* 2.00000 *)

Derivative[0, 1][StieltjesGamma][0, a]
(* 0.324 *)

d // FullForm
(* 0.32399522609896337580027385456880978489`3.339102855094484 *)

Precision[d]
(* 3.3391 *)

Here, one would expect that at least 0.32 are correct digits. But in fact, the exact result is 1 - Pi^2/6 that is -0.644934.... No correct digits, even the sign is wrong.

$\endgroup$
6
  • 3
    $\begingroup$ The authors of the book The SIAM 100-Digit Challenge gave Mathematica (v5) solutions to the 10 problems as well as other implementations. Code here. $\endgroup$
    – Michael E2
    Commented Jun 24, 2013 at 16:09
  • 2
    $\begingroup$ You'll probably find Oleksandr Pavlyks screencast on Mathematical Numerics and Special Functions interesting. $\endgroup$
    – ssch
    Commented Jun 24, 2013 at 16:31
  • 1
    $\begingroup$ +1 toward your first "Good question" badge. $\endgroup$
    – Mr.Wizard
    Commented Jun 24, 2013 at 20:14
  • 1
    $\begingroup$ BTW, I just realized I could write StieltjesGamma[0, #]&'[a] instead of a more verbose Derivative[0, 1][StieltjesGamma][0, a]. $\endgroup$ Commented Jun 24, 2013 at 23:36
  • 4
    $\begingroup$ Probably, you should never expect that NIntegrate returns only provably correct digits. $\endgroup$
    – TauMu
    Commented Jun 27, 2013 at 7:09

4 Answers 4

69
+100
$\begingroup$

Control the Precision and Accuracy of Numerical Results

This is an excellent question.

Of course everyone could claim highest accuracy for their product.

To deal with this situation there exist benchmarks to test for accuracy.

One such benchmark is from NIST. This specific benchmark deals with the accuracy of statistical software for instance.

The NIST StRD benchmark provides reference datasets with certified computational results that enable the objective evaluation of statistical software.

In an old issue of The Mathematica Journal, Marc Nerlove writes elaborately about performing the linear and nonlinear regressions using the NIST StRD benchmark (and Kernel developer Darren Glosemeyer from WRI discussing results using Mathematica version 5.1).


Numerically unstable functions:

But this is only one part of story. OK. There exist benchmark for statistical software etc., but what happens if we take some functions that are numerically unstable?

Stan Wagon has several examples of inaccuracies and how to deal with them in his book Mathematica in Action, which I can only warmly suggest. I have it now for (the latest edition) several years and everytime there is something new to discover with Mr. Wagon.

Let's take, for instance a numerical unstable Maclaurin polynomial of $\sin x$:

poly = Normal[Series[Sin[x], {x, 0, 200}]];
Plot[poly, {x, 0, 100}, PlotRange -> {-2, 2}, 
    PlotStyle -> {Thickness[0.0010], Black}]

The result this we can see that the result breaks down at ~40:

enter image description here

If we take one value x = 60 and perform a division we get a result back:

N[poly /. x -> 60] ==> -0.304811

Inserting the approximate real number 60.; there occurs a roundoff error:

poly /. x -> 60. ==> -4.01357*10^9

But inserting the number 60 (without the period); there is no problem at all:

ply /. x -> 60 ==> -((3529536438455<<209>>9107277890060)/(1157944045943<<210>>4588491415899))

The use of machine precision (caused by the decimal point) leads to an error:

10^17 + 1./100 - 10^17 ==> 0.

Machine precision is $53 \log_{10}(2) = 15.9546$.

This is the exact moment where N comes into play. We have to increase the precision:

poly /. x -> N[60,20] ==> 0. x 10^7

Still not good enough, because this number has no precision at all. So, let's increase the precision again:

poly /. x -> N[60,200] ==> -0.9524129804151562926894023114775409691611879636573830381666715331536022870514582375567159979758451142049758239018693823215314740415313661058559273332324475257579234995809519

This looks much better. If we impose the precision in our prior plot:

Plot[poly, {x, 0, 100}, PlotRange -> {-2, 2}, 
    PlotStyle -> {Thickness[0.0010], Black}, WorkingPrecision -> 200]

enter image description here

Not ideal, since in order to get an accurate result, we need to know what precision we need. There are numerical results which tend to lose precision during several iterations. Luckily there is some salvation in form of the Lyapunov exponent (denoted $\lambda$), which can quantify the loss of precision.

Conclusion:

What I've learned from here is, that it is a bad idea to mix small numbers with big ones in a machine precision environment. This is where Mathematica's adaptive precision comes into play.


Mathematica precision handling

Let's investigate further about precision handling inside Mathematica.

If we want to calculate $\sin(10^{30})$ in Mathematica we get:

N[Sin[10^30]] ==> 0.00933147

Using WolframAlpha we get:

WolframAlpha["Sine(10^30", {{"DecimalApproximation", 1}, "Content"}] ==>  - 0.09011690191213805803038642895298733027439633299304...

The result we get from our numerical workhorse is simply the wrong answer and this is getting worse if we increase the exponent.

(The guys at WolframAlpha seem to do it somewhat differently...but what?)

If we take $10^{30}$ and put turn this into a software real with $MachinePrecision as the actual precision we get 0 as the result, with the precision 0. This result is useless. Luckily we do know that it is indeed.

Here the adaptive precision comes into play.

The adaptive precision is controlled through the system variable $MaxExtraPrecision (default value is 50).

Let's say we want to compute $\sin(10^{30})$ but with a precision of 20 digits:

N[Sin[10^30], 20] ==> -0.090116901912138058030

Ah! We're getting close to the WolframAlpha engine!

If we ask for $\sin(10^{60})$ the result is:

N[Sin[10^60], 20] ==> N::meprec: Internal precision limit
                      $MaxExtraPrecision = 50.` reached while evaluation
     Sin[1000000000000000000000000000000000000000000000000000000000000]. >>
     Out[105]= 0.8303897652

We run into problems, since the adaptive algorithm only adds 50 digits for extra precsion. But, luckily, the extra precision is controlled through $MaxExtraPrecision, which we're allowed to change:

$MaxExtraPrecision = 200; N[Sin[10^60], 20] ==> 0.83038976521934266466

Addendum (Michael E2):

Note that N[Sin[10^30]] does all the computation in MachinePrecision without keeping track of precision; however N[Sin[10^30], n] does keep track and will give an accurate answer to precision n. (WolframAlpha probably uses something like n = 50.) Also specifying the precision of the input to be, say, 100 digits,N[Sin[10^60`100], 20] will use 100-digit precision calculations internally and return the same answer as above to 20 digits of precision, provided as in this case 100 digits is enough to give 20. (Added at the request of @stefan.)

Conclusion

Equipped with that knowledge we could define functions that use adaptive precision to get an accurate result.


Precision and accuracy

It is not that Mathematica loses precision, but in your definition of a you'll lose precision in the first place.

Let's first talk about precision and accuracy.

Basically the mathematical definition of precision and accuracy is as follows:

Suppose representation of a number $x$ has an error of size $\epsilon$. Then the accuracy of $x \pm \epsilon/2$ is defined to be $-\log_{10}|\epsilon|$ and its precision $-\log_{10}|\epsilon/x|$.

With these definitions we can say that a number $z$ with accuracy $a$ and precision $p$ will lie with certainty in the interval:

$$\left(x-\frac{10^{-a}}{2},\frac{10^{-a}}{2}\right)=\left(x-\frac{10^{-p} x}{2},\frac{10^{-p} x}{2}+x\right)$$

According to these definitions the following relation holds between precision and accuracy:

$\operatorname{precision}(x)=\operatorname{accuracy}(x)+\log_{10}(|x|),$

where the latter is called the scale of the number $x$.

We can check if this identity holds:

Function[x, {Precision[x], Accuracy[x] + Log[10, Abs[x]]}] /@
         {N[1, 100], N[10^100, 30]}

==> {{100.,100.},{30.,30.}} (* qed *)

Let's define a function for both precision and accuracy:

PA[x_] := {Precision[x], Accuracy[x]}

Now let's look at your definition of a:

a = 1`7

PA[a] ==> {7., 7.}

d = Derivative[0, 1][StieltjesGamma][0, a] ==> -1.6450

PA[d] ==> {5.15586, 4.93969}

You've lost precision!

You defined a to have a precision and an accuracy of 7.

But what is the precision and accuracy if you turn a into a symbol using machine precision:

a = 1.

PA[a] ==> {MachinePrecision, 15.9546}

This is a gain in precision obviously. Now let's call your canonical examples:

d = Derivative[0, 1][StieltjesGamma][0, a]

==> -1.64493

Which is the exact result of $-\frac{\pi ^2}{6}$.

The precision and accuracy of d is:

PA[d] ==> {MachinePrecision, 15.7384}

Perfect.

Now let's redefine your a to be 2. instead of 2`6:

a = 2.

PA[a] ==> {MachinePrecision, 15.6536}

d = Derivative[0, 1][StieltjesGamma][0, a]

==> -0.644934

Which is the exact result of $1 - \frac{\pi ^2}{6}$

PA[d] ==> {MachinePrecision, 16.1451}

Conclusion

Dealing with numerical computing is dealing with loss of precision. It seems that Mathematica varies the Precision depending on the numerical operation being performed and the Precisions are more pessimistic than optimistic, which is actually quite good.

In most calculations, one typically loses precision, but with an appropriate starting value you can gain precision as well.

The general rule for the usage of high-precision numbers is:

If you want to gain high precision you need to use high-precision numbers in your expression to be calculated. Consequently, every time you need a high-precision result you must take care that the starting expression has sufficient precision.

There exists an exception to the above rule. If you use machine-precision arithmetic in expressions and the numbers are getting bigger than $MaxMachineNumber, Mathematica will switch automatically to high-precision numbers. If this is the case the rules apply as described in my Edit 2.

P.S.:

This was one of the questions I really like, since I know now more about that topic than before. Maybe one of the WRI/SE jedi can join the party to provide even more insights on this, possibly more than I would ever be able to provide.

$\endgroup$
6
  • $\begingroup$ @MichaelE2 Thank you for your helpful input. I'd like to invite you, if you can afford the time, to update the post and enter your findings with your attribution. If this is ok for you. Or post it in a new post, if this is to circuitous for you. But thanks again for your input. $\endgroup$
    – Stefan
    Commented Jun 24, 2013 at 15:14
  • $\begingroup$ Thanks. I added it, with correct formatting :), as an add-on to "Edit 2". You may wish to roll back, or re-edit it yourself, if you would like to incorporate it more effectively. That would be fine with me. (I'll delete my initial comment soon, I think.) $\endgroup$
    – Michael E2
    Commented Jun 24, 2013 at 15:54
  • $\begingroup$ @MichaelE2 nice job :) i appreciate your input indeed. $\endgroup$
    – Stefan
    Commented Jun 24, 2013 at 19:42
  • 1
    $\begingroup$ @Stefan Thank you for your deep analysis of the subject! The bounty is yours. $\endgroup$ Commented Jul 2, 2013 at 18:32
  • $\begingroup$ @VladimirReshetnikov thank you for accepting it :) I was about to write an addendum on NumericalMath`$NumberBits function. This function shows how Mathematica simulates interval arithmetic by constantly maintaining a few more digits than needed... $\endgroup$
    – Stefan
    Commented Jul 2, 2013 at 18:36
5
$\begingroup$

This isn't an answer (yet) but it was too long for a comment. Here is an extended example where the quoted precision does not appear to be true:

f1 = Derivative[0, 6][StieltjesGamma][0, #] &;
f2 = {Accuracy@#, Precision@#, InputForm@#} &;

f1 @ 1`15
f2 @ %

N[f1 @ 1, 10]
f2 @ %
725.59

{2.29961, 5.1603, 725.59278417719148802796650335754143851756`5.1603018938964}

726.0114797

{7.13906, 10., 726.01147971477215163394242248090753970333`10.}
$\endgroup$
5
$\begingroup$

Update

As @MichaelE2 points out, perhaps I didn't directly answer your question. I think the short answer is that you can rely on precision of the output if you are using adaptive precision control, and you can't rely on the precision of the output if you are not using adaptive precision control. In your examples you are not using adaptive precision control, so your results are not surprising.

Extended response

If you have a function, and you insert an inexact number, then Mathematica has very little control over the precision/accuracy of the output. Here is an extremely simple example:

Precision[1.0001`5 - 1]

0.999957

You started off with a precision of 5, and then Mathematica returned a result with a precision of 1. The loss of precision is due to subtractive cancellation. If instead you insert an exact number and numericize, then Mathematica can use adaptive precision control to insure that the output has the precision requested during the numericization process. Adaptive precision control basically means that while numericizing an input, exact numbers can be numericized to much higher precision when needed, controlled by the global variable $MaxExtraPrecision. More schematically:

Precision[f[N[exact, prec]]] ->  ?
Precision[N[f[exact], prec]] -> prec

So, in your question, you are doing the former, and hence you are circumventing Mathematica's adaptive precision control, and so, just like with MachinePrecision computations, the result may not have any correct digits. Here's another simple example. Suppose your function is:

f[x_] := 1/1000 - Sin[x];

If we feed f an inexact number close to 1/1000, then due to subtractive cancellation, the output will have much less precision than the input:

Precision @ f[N[1/1000, 10]]

3.22185

On the other hand, if you feed an exact number, and then numericize, the evaluation will use adaptive precision control:

Precision @ N[f[1/1000], 10]

10.

We can use the following trick to see adaptive precision control in action:

Clear[f]
f[x_Real] := 1/1000 - Sin[x]

TracePrint[
    N[f[1/1000], 10],
    _f,
    TraceInternal->True,
    TraceAction->(Print[InputForm[#]]&)
]
Precision[%]

HoldForm[f[1/1000]]

HoldForm[f[1/1000]]

HoldForm[f[0.001`10.]]

HoldForm[f[0.001`29.93398985691138]]

1.666666583*10^-10

10.

You can see that N first feeds f an inexact number with precision 10, and as we saw earlier, this is not sufficient to produce an answer with the requested precision. Then, N feeds f an inexact number with precision 30, and this was sufficient.

The moral of the story is N[f[exact], prec] should produce a result with precision prec, while f[N[exact, prec]] need not.

$\endgroup$
1
  • $\begingroup$ Concerning calculations like your first example, isn't the OP's question, when can we be sure the first digit is correct if the Precision is at least 1? In the OP's examples, the Precision is greater than the actual precision of the result. My read of the OP is that the question is about the accuracy of Precision, not about how to achieve a certain precision. $\endgroup$
    – Michael E2
    Commented Feb 8, 2018 at 18:22
1
$\begingroup$

This is really a big topic at Wolfram Inc. As this post states compare-their-asymptotic-behavior-using-the-order-notation. Wolfram Language did improve much at the methodologies of making deviation of values more present to the users. There are many questions and answers in this community about big O notation and alike. Series had the big O notation very early through the sequels of Mathematica and Wolfram Language.

Since version 12 there is a bunch of new built-ins starting with the Asymptotic pre word like AsymptoticSolve doing mathematics at critical points like zero, singularities of polynomial and nonpolynomial order. So the range where mathematical correctness and solution ability exist and can be computed with Mathematica is extended to what is thought at present time universities.

As computational mathematicians accept today the approximate is accepted as exactness of its own value. So there is much more than function approximation or calculation of values at chosen points.

An example is PadeApproximation which uses rational functions as approximations to functions that are inaccessible to be approximated with plain polynomials. There are even built-ins like GeneralMiniMaxApproximation that allow setting a brake to the validity domain of functional relations. This is a range of approximation theory making more clear, why the Taylor expansion has to fail on transcendent functions because this behavior of very much located deviation from exactness is much more present.

It is used despite this singularity of unavoidable order to extend exactness ranges in domains.

Since the introduction of big o notation the users of the Wolfram Language can compute for themself with exactness measures beyond for example. Beyond for example Accuracy. The definitions of Accuracy are

... gives the effective number of digits to the right of the decimal point in the number x.

... gives a measure of the absolute uncertainty in the value of x.

For exact numbers such as integers, Accuracy is Infinity.

For any approximate number x, Accuracy is equal to Precision - RealExponent. The core set of built-ins for correctness in Wolfram Language is defined in RepresentationOfNumbers. So the logic is in the section testing for types. And the access to the internal is in Internal Representation.

In this, Your posed question targets the definition of [$MaxPrecision][11] at most. The default value of $MaxPrecision is Infinity. Which is not true for most calculations, and computations with Wolfram Language. But this is a general aim that is possible. The issue given for this parameter is: When $MaxPrecision is set, some computations may not get correct results! This is meant a priori so before the operation is executed.

It might be possible to get past the limit but needs skills that are not in general implemented in Mathematica or Wolfram Language. There are many pages on Mathematica at Stackexchange addressing real advances. As the tutorial ArbitraryPrecisionNumbers shows not in all cases is the extension of exactness for the input better than working with inexact input. This is where consideration have to start. Conduction of the numerical experiment, if one is not sure about exactness, is the next step.

The guide NumbersWithUncertainty starts with the foundation of error propagation. The guide PrecisionAndAccuracyControl continues this into the control of output for accuracy and precision. These two terms substitute or extent the concept of correct or correctness far. For physicists, this can be done much further than experimental precision and exactness since numerical mathematics is targeted toward infinite correctness.

To return to a professional path the TimeMeasurementAndOptimization may be used to connect the desire, the need for correctness to business style payable toughness. Not all that can be computed is in need or necessary beyond the investment of time. This is for many people incorporated in PerformanceGoal for what is really reaching my customers as value or an increase in value.

So SystemOptions is the next place where to do settings for correctness. This page open another perspective on exactness: LowLevelSystemSpelunking. The original wording is:

"The Wolfram System is a large and complex software system. Although strongly not supported for production purposes, it is sometimes instructive to "spelunk" in the system, looking at the low-level organization and obscure features." So searching the Mathematica documentation for obscure leads to the image built-in Blur. This is used to hide details on images and indeed one can find examples where obscuring is useful. The example of choice is the galloping horse. Where the stroboscopical pictures of a galloping horse are overlayed to prove in one view that there are phases in the gallop of a horse where the horse has no contact to the ground.

First@Image3DSlices[Blur[Image3D@anim, {3, 0, 0}], {12}]

gallop study of a horse with a person on it

One glance and everything in question is shown. But the picture is highly processed and far from reality and therefore obscure. This did make a British man famous and rich and nowadays it is common knowledge but the representation is only simple in the obscured picture. And so the horse flies in the gallop and it is never considered a flying species. And that is what mathematical correctness is too.

$\endgroup$

Your Answer

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