## Wednesday, June 11, 2008

### Basic implementation of significance arithmetic

This .py file contains a work-in-progress implementation of significance arithmetic, using mpmath as a base. The main class is currently called "num" (sorry, I just haven't come up with a good name yet). An instance can be created from a float or int value. A second number that specifies the accuracy (measured in decimal digits) can be passed; floats are assumed by default to have an accuracy of 53 bits or 15.95 decimal digits. The method .accuracy(b) gives the estimated number of accurate digits in the given base (default b = 10).

`>>> num(1.2345678901234567890123456789)1.23456789012346>>> _.accuracy()15.954589770191001>>> num(1.2345678901234567890123456789, 5)1.2346>>> _.accuracy()5.0>>> num(1.2345678901234567890123456789, 25)1.234567890123456690432135>>> _.accuracy()25.0`

In the last example, the fact that the input float has a limited accuracy as a representation of the entered decimal literal becomes visible. (Support for passing an exact string as input to avoid this problem will be added.)

The accuracy is a measure of relative accuracy (or relative error). A num instance with value u and accuracy a represents the interval u · (1 + ξ 2-a) for values of ξ between -1 and 1. The relative error, in the traditional numerical analysis meaning of the term, is given by 2-a and the absolute error is given by |u| · 2-a. In other words, the accuracy is a logarithmic (base 2) measure of error; in some cases, it is more natural to consider the (logarithmic) error, given by -a.

It can usually be assumed that a > 0. When this property does not hold, not even the sign of the number it is supposed to represent is known (unless additional assumptions are made). Such extremely inaccurate numbers will need special treatment in various places; when they result from higher operations (such as when the user asks for N digits of an expression in SymPy), they should probably signal errors.

How do errors propagate? The simplest case is that of multiplication. If x = u · (1+e) and y = v · (1+f), then x · y = u · v · (1+e) · (1+f) = (u · v) · (1 + e + f + e·f). Therefore the cumulative error is given by e + f + e·f. Expressed in terms of the logarithmic representation of the errors, 2-a = e and 2-b = f, the final accuracy is given by c = −log2(2-a + 2-b + 2-a-b). This expression can usually be approximated as c = min(a,b), which is the rule I have implemented so far. (Further testing will show whether it needs to be refined).

This analysis does not account for the error due to roundoff in the floating-point product u · v. Such errors can be accounted for either by adding an extra term to the error sum or by setting the arithmetic precision slightly higher than the estimated accuracy. I'm currently taking the latter approach.

Division is similar to multiplication.

Addition and subtraction is a bit harder, as these operations need to translate the relative errors of all the terms into a combined absolute error, and then translate that absolute error back to a relative accuracy (of the final sum). It is important to note that the accuracy of a sum can be much greater or much lesser than the minimum accuracy of all the terms.

Translating between relative and absolute error associated with a number involves knowing its magnitude, or rather the base-2 logarithm thereof. This is very easy with mpmath numbers, which are represented as tuples x = (sign, man, exp, bc) where bc is the bit size of the mantissa. The exact magnitude of x is given by log2 x = log2 man + exp, and this quantity can be approximated closely as exp+bc (although doing so admittedly does not pay off much since a call to math.log is inexpensive in Python terms).

Perhaps the most important type of error that significance arithmetic catches is loss of significance due to cancellation of terms with like magnitude but different sign. For example, 355/113 is an excellent approximation of pi, accurate to more than one part in a million. Subtracting the numbers from each other, with each given an initial accuracy of ~15 digits, leaves a number with only 9 accurate digits:

`>>> pi = num(3.1415926535897932)>>> 355./113 - pi2.66764189e-7>>> _.accuracy()9.0308998699194341`

Compare this with what happens when adding the terms:

`>>> 355./113 + pi6.28318557394378>>> _.accuracy()15.954589770191001`

Adding an inaccurate number to an accurate number does not greatly reduce accuracy if the less accurate number has small magnitude:

`>>> pi + num(1e-12,2)3.1415926535908>>> _.accuracy()14.342229822223226`

Total cancellation currently raises an exception. Significance arithmetic requires special treatment of numbers with value 0, because the relative error of 0 is meaningless; instead some form of absolute error has to be used.

I have also implemented real-valued exponentiation (x**y, exp(x) and sqrt(x)). These operations are currently a little buggy, but the basics work.

Exponentiation reduces accuracy proportionally to the magnitude of the exponent. For example, exponentiation by 1010 removes 10 digits of accuracy:

`>>> from devel.num import *>>> pi = num(3.1415926535897932)>>> pi.accuracy()15.954589770191001>>> pi ** (10**10)8.7365e+4971498726>>> _.accuracy()5.954589770191002`

Contrast this with the behavior of mpmath. With the working precision set to 15 digits, it will print 15 (plus two or three guard) digits of mantissa. Re-computing at a higher precision however verifies that only the first five digits were correct:

`>>> from mpmath import mp, pi>>> mp.dps = 15>>> pi ** (10**10)mpf('8.7365179634758897e+4971498726')>>> mp.dps = 35>>> pi ** (10**10)mpf('8.7365213691213779435688568850099288527e+4971498726')`

One might wonder where all the information in the input operand has gone. Powering is of course implemented using binary expontiation, so the rounding errors from repeated multiplications are insignificant.

The answer is that exponentiation transfers information between mantissa and exponent (contrast with a single floating-point multiplication, which works by separately multiplying mantissas and adding exponents). So to speak, exponentiation by 10n moves n digits of the mantissa into the exponent and then fiddles a little with whatever remains of the mantissa. Logarithms do the reverse, moving information from the exponent to the mantissa.

This is not something you have to think of often in ordinary floating-point arithmetic, because the exponent of a number is limited to a few bits and anything larger is an overflow. But when using mpmath numbers, exponents are arbitrary-precision integers, treated as exact. If you compute pi raised to 1 googol, you get:

`>>> num(3.1415926535897932) ** (10**100)1.0e+4971498726941338374217231245523794013693982835177254302984571553024909360917306608664397050705368816>>> _.accuracy()-84.045410229808965`

The exponent, although printed as if accurate to 100 digits, is only accurate to 15. Although the accuracy is reported as negative, the number does have a positive "logarithmic accuracy". So in contexts where extremely large numbers are used, some extra care is needed.

A counterintuitive property of arithmetic, that a direct implementation of precision-tracking floating-point arithmetic fails to capture, is that some operations increase accuracy. Significance arithmetic can recognize these cases and automatically adjust the precision to ensure that no information is lost. For example, although each input operand is accurate to only 15 decimal places, the result of the following operation is accurate to 65:

`>>> num(2.0) ** num(1e-50)1.0000000000000000000000000000000000000000000000000069314718055995>>> _.accuracy()65.783713593702743`

This permits one to do things like

`>>> pi = num(3.1415926535897932)>>> (pi / (exp(pi/10**50) - 1)) / 10**501.0>>> _.accuracy()15.65355977452702`

which in ordinary FP or interval arithmetic have a tendency to cause divisions by zero or catastrophic loss of significance, unless the precision is manually set high enough (here 65 digits) from the start.

Automatically-increasing precision is of course a bit dangerous since a calculation can become unpexpectedly slow in case the precision is increased to a level much higher than will be needed subsequently (e.g. when computing 15 digits of exp(10-100000) rather than exp(-100000)-1). This feature therefore needs to be combined with some user-settable precision limit.

I'll finish this post with a neat example of why significance or interval arithmetic is important for reliable numerical evaluation. The example is due to Siegfried M. Rump and discussed further in the paper by Sofroniou and Spaletta mentioned in an earlier post. The problem is to evaluate the following function for x = 77617 and y = 33096:

`def f(x,y):    return 1335*y**6/4 + x**2*(11*x**2*y**2-y**6-121*y**4-2) + \        11*y**8/2 + x/(2*y)`

Directly with mpmath, we get (at various levels of precision):

`>>> from mpmath import mp, mpf, nstr>>> for i in range(2,51):...     mp.dps = i...     print i, nstr(f(mpf(77617),mpf(33096)),50)...2 1.1718753 2596148429267413814265248164610048.04 1.1726074218755 1.1726036071777343756 1.17260384559631347656257 1.17260393500328063964843758 -9903520314283042199192993792.09 -1237940039285380274899124224.010 -154742504910672534362390528.011 9671406556917033397649408.012 1.1726039400532499712426215410232543945312513 75557863725914323419136.014 -9444732965739290427392.015 -1180591620717411303424.016 1.172603940053178639413289374715532176196575164794917 -18446744073709551616.018 1152921504606846977.2519 144115188075855873.17187520 1.172603940053178631859449551101681752385275103733921 1125899906842625.172603845596313476562522 1.172603940053178631858840746170842724016569746936523 1.17260394005317863185883412872594229979517077566724 1.172603940053178631858834955906554852822845647075725 1.172603940053178631858834904207766568258615967612626 8589934593.172603940053178625535501566901075420901227 1073741825.172603940053178631823874167317001138144428 1.172603940053178631858834904510689155863484500890729 1.172603940053178631858834904520155486726136642555730 1.172603940053178631858834904520155486726136642555731 1.172603940053178631858834904520180138629424799174632 1.17260394005317863185883490452018322011733581875233 1.172603940053178631858834904520183797896319134922734 1.172603940053178631858834904520183701599821915560935 1.17260394005317863185883490452018370761835299177136 -0.8273960599468213681411650954798162920054888159658437 -0.8273960599468213681411650954798162920054888159658438 -0.8273960599468213681411650954798162919996113442117339 -0.8273960599468213681411650954798162919990603312347840 -0.8273960599468213681411650954798162919990373723607441 -0.8273960599468213681411650954798162919990330675718642 -0.8273960599468213681411650954798162919990330675718643 -0.8273960599468213681411650954798162919990331124134144 -0.8273960599468213681411650954798162919990331152160145 -0.8273960599468213681411650954798162919990331157414946 -0.8273960599468213681411650954798162919990331157852847 -0.8273960599468213681411650954798162919990331157852848 -0.8273960599468213681411650954798162919990331157844349 -0.8273960599468213681411650954798162919990331157843950 -0.82739605994682136814116509547981629199903311578438`

Remarkably, the evaluations are not only wildly inaccurate at low precision; up to 35 digits, they seem to be converging to the value 1.1726..., which is wrong!

Significance arithmetic saves the day:

`>>> for i in range(2,51):...     try:...         r = f(num(77617,i),num(33096,i))...         s = str(r)...         a = r.accuracy()...         print i, str(s), a...     except (NotImplementedError):...         continue...2 2.0e+33 -2.515449934965 6.0e+29 -2.826779887266 -8.0e+28 -2.729869874267 -1.0e+28 -2.632959861258 -1.0e+27 -2.536049848249 -2.0e+26 -2.4391398352310 1.0e+25 -2.6432598178912 2.0e+23 -2.4494397918713 -9.0e+21 -2.6535597745314 -1.0e+21 -2.5566497615216 -2.0e+19 -2.362829735517 1.0e+18 -2.5669497181618 -7.0e+16 -2.7710697008120 1.0e+15 -2.577249674821 -1.0e+14 -2.4803396617925 9.0e+9 -2.6947596010926 1.0e+9 -2.5978495880827 -1.0e+8 -2.5009395750735 -8.0e-1 -2.9297794536636 -8.0e-1 -1.9297794536637 -8.0e-1 -0.92977945366238 -8.0e-1 0.070220546338439 -8.0e-1 1.0702205463440 -8.3e-1 2.0702205463441 -8.27e-1 3.0702205463442 -8.274e-1 4.0702205463443 -8.274e-1 5.0702205463444 -0.827396 6.0702205463445 -0.8273961 7.0702205463446 -0.82739606 8.0702205463447 -0.82739606 9.0702205463448 -0.8273960599 10.070220546349 -0.82739605995 11.070220546350 -0.827396059947 12.0702205463`

I had to wrap the code in try/except-clauses due to num(0) not yet being implemented (at many of the low-precision stages, but fortunately for this demonstration not all of them, there is not unpexpectedly complete cancellation).

The accuracy is reported as negative until the arguments to the function are specified as accurate to over 38 digits, and at that point we see that the printed digits are indeed correct. (The values do not exactly match those of mpmath due to slightly different use of guard digits.)

Interestingly, the num class seems to greatly overestimate the error in this case, and I'm not yet sure if that's inherent to evaluating Rump's particular function with significance arithmetic or due to implementation flaws. It's of course better to overestimate errors than to underestimate them.

Numerical evaluation of a SymPy expression can be seen as converting the expression to a function of its constants (integers and special numbers like pi) and evaluating that function with the constants replaced by significance floating-point numbers. In practice, I will probably implement numerical evaluation using a recursive expression tree walker. This way forward error analysis can be used to efficiently obtain the user-requested level of accuracy at the top of the expression. Subevaluations can be repeated with increased precision until they become accurate enough; precision only needs to be increased for parts of an expression where it actually matters.

My next goal, besides polishing the existing features (and implementing num(0)), is to implement more functions (log, trigonometric functions), comparison operators, and writing some tests. With trigonometric functions, all the arithmetic operations for complex numbers will be straightforward to implement.

Anonymous said...

http://www2.hursley.ibm.com/decimal/decifaq4.html#signif

(and, in particular, Delury's paper referred to there).

Mike Cowlishaw

Fredrik Johansson said...

Thanks for the input, Mike.

Significance arithemtic is certainly bad in some situations. A rather serious problem that I think your FAQ neglects to mention is that it is useless for "self-correcting" numerical algorithms such as Newton's method.

What I'm trying to accomplish right now is mainly reliable conversion of exact mathematical expressions into accurate numerical values, similar to what Mathematica does.

For working with measurements in a more rigorous way, I would like to some day extend SymPy's statistics module so that you can do arithmetic with (arbitrary) statistical distributions.

Anonymous said...

OK, thanks -- good point re. converging algorithms.

Mike

newBBie said...

I'm trying to find the correct formula for "Rump's Example" (and free of typo errors). Please see slide 6 of this: https://hpcrd.lbl.gov/SciDAC08/files/presentations/ARendell.SciDAC.pdf

It doesn't agree with your formula for Rump, and I believe the slide has other problems.

Thanks!