Interval Arithmetic

1. Interval Arithmetic Kwang Hee Ko School of Mechatronics Gwangju Institute of Science and Technology

2. Interval Numbers • An interval number is defined as • An ordered pair of real numbers denoted by [a,b] • It is a set of real numbers x such that • [a,b] = {x|a ≤ x ≤ b} • Operations in set theory can be applied. • Two intervals [a,b] and [c,d] are said to be equal if • a=c and b=d • The intersection of two intervals is empty of [a,b]∩[c,d] = 0 if either • a > d or c > b • Otherwise [a,b]∩[c,d]=[max(a,c),min(b,d)] • The union of the two intersecting intervals is • [a,b]∪[c,d]=[min(a,c),max(b,d)]

3. Interval Numbers • An order of intervals is defined by • [a,b]<[c,d] if and only if b < c • Width: w([a,b]) = b-a • Magnitute: |[a,b]| = max(|a|,|b|) • Examples • [2,4]∩[3,5] = [3,4] • [2,4]U[3,5] = [2,5] • |[-7,-2]| = 7

4. Interval Arithmetic • Exact-Interval Arithmetic: an extension of real arithmetic. We assume that endpoints are computed with infinite precision. • Basic arithmetic operations: +,-,·,/ [a,b]*[c,d] = {x*y| a ≤ x ≤ b, c ≤ y ≤ d} 0 ∈ [c,d]

5. Interval Arithmetic • Examples • [2,4]+[3,5] = [2+3,4+5] = [5,9] • [2,4]-[3,5] = [2-5,4-3] = [-3,1] • [2,4]*[3,5]=[min(2*3,2*5,4*3,4*5),max(2*3,2*5,4*3,4*5)]=[6,20] • [2,4]/[3,5]=[min(2/3,2/5,4/3,4/5),max(2/3,2/5,4/3,4/5)]=[2/5,4/3]

6. Interval Arithmetic • Interval arithmetic is associative, commutative. • NOT distributive but subdistributive [a,b]ᆞ[e,f]

7. Interval Arithmetic • Example of being subdistributive • [1,2]*([1,2]-[1,2]) = [1,2]*([-1,1])=[-2,2] • [1,2]*[1,2] – [1,2]*[1,2] = [1,4]-[1,4] = [-3,3] • [-2,2] ⊆ [-3,3]

8. Interval Arithmetic • Inclusion Monotonic • If I ⊂ K, J ⊂ L • I + J ⊂ K + L • I – J ⊂ K – L • IJ ⊂ KL • I/J ⊂ K/L (if 0 ∈ L)

9. Rounded Interval Arithmetic • If floating point arithmetic is used to evaluate these interval arithmetic equations, there is no guarantee that the roundings of the bounds are performed conservatively. • Rounded interval arithmetic (RIA) ensures that the computed end points always contain the exact interval.

10. IEEE 754 Standard format of Floating Point Arithmetic • Units-in-last-place, ulpl, ulpu for each separate floating point number. • Concept • To calculate ulp it is necessary to extract the integer value of the exponent from the binary representation. • Recall that the value of the significand B of a double precision number X is • X = (-1)s2EB with B = 1+b12-1+b22-2+…+b522-52 • The value of the least significant bit b52 is 2-52. • Thus, the value of ulp is 2E2-52 = 2E-52.

11. IEEE 754 Standard format of Floating Point Arithmetic • Using the C standard functions frexp() and ldexp(). • Because of the use of standard library functions, this implementation is slow.

12. IEEE 754 Standard format of Floating Point Arithmetic • To avoid using the library functions and to construct the ulp directly… • The biased exponent e occupies bits 52 through 62. • If we could manipulate the binary representation as a 64-bit integer, we could extract e by dividing by 252. • Right-shift the bit pattern by 52 bits, placing e in bits 10 through 0. • The sign bit would then occupy bit 11, could be removed by performing a bitwise logical AND with the 64-bit mask 0…011111111111. • Most commercial processors and programming languages do not support 64-bit integers; generally, only 8, 16 and 32-bit integers are available. • To overcome this limitation, we can overlay the storage location of the 640bit double precision number with an array of four 16-bit integers (short type)

13. IEEE 754 Standard format of Floating Point Arithmetic • In C or C++ this can be accomplished using the union data structure • typedef union { double dp; /* the 64-bit double precision value */ unsigned short sh[4] /* overlay an array of 4 16-bit integers*/ } Double • After the assignment • Double x; • Double D; • D.dp = x;

14. IEEE 754 Standard format of Floating Point Arithmetic • In C or C++ this can be accomplished using the union data structure • The exponent of the variable x can be extracted from D.sh[0], whose 16 bits contain • the sign bit s (bit 15), • the 11-bit biased exponent e (bits 14 through 4) • The 4 most significant bits b1b2b3b4 of the mantissa m (bits 3 through 0) • The biased exponent e can be extracted from D.sh[0] by performing a bitwise logical AND with the 16-bit mask 0111111111110000 to zero-out the sign bit and the four most significant bits of the mantissa and then right-shifting e by 4 bits. Then the unbiased exponent E = e-1023.

15. IEEE 754 Standard format of Floating Point Arithmetic • Algorithm

16. IEEE 754 Standard format of Floating Point Arithmetic • Comparison of two different implementations • The timings (in CPU seconds) were taken on a 100MHz RISC processor. • The values are the accumulated CPU times to perform 100,000 calculations of the ulp of various representative values of X.

17. IEEE 754 Standard format of Floating Point Arithmetic • Software/hardware roundings in implementation of rounded interval arithmetic • The software rounding method is computationally more expensive than hardware rounding, requiring an extra addition and subtraction and the computation of the ulp of two values. • The software rounding method extends the upper and lower bounds of the interval during every arithmetic operation. • Hardware rounding is superior to software rounding. • Set a CPU flag for appropriate roundings. • The hardware rounding method only extends the bounds when the result of the operation cannot be exactly represented, producing tighter interval bounds.

18. IEEE 754 Standard format of Floating Point Arithmetic • Whenever an exact value cannot be represented in a given precision, rounding occurs. • IEEE 754 Standard specifies four rounding modes: round to nearest, round to positive infinity, round to negative infinity and round to zero. • As default, the standard sets the round to nearest mode.

19. IEEE 754 Standard format of Floating Point Arithmetic • MS Windows (Visual C++) • Linux (Refer to ‘fpu_control.h’ file) • _FPU_SETCW(arg) • -infinity: arg = 0x077f • +infinity: arg = 0x0b7f • Nearest : arg = …

20. RIA Implementation • Use C++

21. RIA Implementation • We can define various additional operations and functions other than basic arithmetic ones in terms of interval arithmetic. • Sqrt, exp, sin, cos, etc. • Disadvantages • Performance is slower than floating point arithmetic. • It tends to suffer overestimation of intervals during computation.

22. Alternative to Interval Arithmetic • Affine Arithmetic • Conceptually similar to interval arithmetic in that affine arithmetic also maintains a range containing the exact value. • But it works better to control the size of a range.

23. Robustness Issues • Suppose we have a degree four planar Bezier curve whose control points are given by • (-0.5,0.0635), (-0.25,-0.0625), (0,0.0625), (0.25,-0.0625), (0.5,0.0625) • It is equivalent to the explicit curve y=x4 (-0.5 ≤ x ≤ 0.5) • The curve intersects with x-axis tangentially at (x,y) = (0,0)

24. Robustness Issues • We translate the curve by +1 in the y direction and translate it back to the original position by moving it by -1/3 three times. • The curve is generally not the same as the original one in floating point arithmetic. • Ex.) We evaluate the curve at parameter value t=0.5 to get (0,0.00035) instead of (0,0). • This is a serious change of topological structure of the curve. • Intersection -> no intersection • Geometric operations in floating point arithmetic is subjected to such errors.

25. Robustness Issues • Floating Point Arithmetic • Limited precision / errors due to rounding : possibility of missing roots

26. Examples in Real World • Some disasters attributable to bad numerical computing in real life • Patriot Missile Failure: in Dharan, Saudi Arabia, on February 25, 1991. This is ultimately attributable to poor handling of rounding errors • Explosion of the Ariane 5 rocket: Just after lift-off on its maiden voyage off French Guiana, on June 4, 1996, it exploded, which was ultimately the consequence of a simple overflow. • The sinking of the Sleipner A offshore platform: Due to inaccurate finite element analysis, the platform sank down to the ocean floor costing nearly one billion dollars.

27. The Patriot Missile Failure • Inaccurate time computation • The time in tenths of second as measured by the system’s internal clock was multiplied by 1/10 to produce the time in seconds. • The calculation was performed using a 24 bit fixed point register. • The value 1/10 is a non-terminating binary expansion, which is chopped at 24 bits. • 1/10 = 1/24 + 1/25 + 1/28 + … • The Patriot missile stores 24 bits only. • An error, 0.000000095, is introduced. • 100 hour operation -> 0.000000095 * 100 * 60 * 60 * 10 = 0.34sec. This is enough time to miss a Scud since during 0.34 sec, the Scud travels more than half a kilometer.

28. Ariane 5 Rocket Explosion • Insufficient Precision • Ariane 5, which is worth 7 billion dollars crashed in less than a minute after lift-off. • A small computer program is to blame for the crash which tries to stuff a 64-bit number into a 16-bit space. • At the bottom of a sequence of events to disaster, conversion of velocity data in a 64-bit format to a 16-bit format triggered an overflow error, making the core system in the rocket shut down.

29. Robustness in Numerical Computation IRoot Finding Kwanghee Ko School of Mechatronics Gwnagju Institute of Science and Technology

30. Motivation • Why do we need a nonlinear polynomial solver? • Many geometric problems are formulated as systems of nonlinear polynomial equations • Distance function • Intersections • Curvature extrema, etc.

31. Solution Methods • Local Solution Methods • Newton-type method • good initial approximation • No assurance that all roots have been found • Global Solution Methods • Algebraic Type Methods • Homotopy (Continuation) Methods • Subdivision Methods

32. Subdivision Methods • Advantages • Efficient and stable • Easy to implement • Disadvantages • No certainty that each root has been isolated • No explicit information about root multiplicities • Example Algorithm • Projected Polyhedron Algorithm

33. Projected Polyhedron Algorithm • Transform algebraic problems (finding roots) to geometric problems (computing intersections). • Use a powerful geometric property, the convex hullproperty in the algorithm. • Input : a system of polynomial equations in Bernstein form. • Output : intervals (regions) which contain roots. Note : Conversion from power basis to Bernstein basis is unstable. So we need to use Bernstein polynomials from the beginning to formulate a problem.

34. Projected Polyhedron Algorithm(Univariate Polynomial Case) f(t) • Make an transformation such that the range t of a function f(t) is from 0 to 1. • Construct a graph (t,f(t)) 0 t1 t2 1 t

35. Projected Polyhedron Algorithm (Univariate Polynomial Case) Construct the convex hull of the Bezier curve

36. Projected Polyhedron Algorithm (Univariate Polynomial Case) Intersect the convex hull with the parameter axis

37. Projected Polyhedron Algorithm (Univariate Polynomial Case) Discard the regions which do not contain roots after applying the de Casteljau subdivision algorithm Eliminate

38. Projected Polyhedron Algorithm (Multivariate Polynomial Case) • x = (x1,x2,…,xm) independent variables • f1(x)=f2(x)=…=fn(x)=0 • Algorithm • Affine Transformation : 0  xi  1, i=1,…,m • Construct a graph for each fj, j=1,…,n. • Project the control polygon of each fj onto m different coordinate (2D) plane, i.e. (x1,xm+1)-plane, …, (xm,xm+1)-plane. • Construct n two dimensional convex hulls in each plane. • Intersect each convex hull with the horizontal axis : x1,…,xm. • Compute intersection of all the intervals obtained at step 5. • If empty, no solution. If non-empty, discard the portions that do not contain roots. • Repeat.

39. Robustness Issues in Root Finding • Floating Point Arithmetic • Limited precision / errors due to rounding : possibility of missing roots

40. Robustness Issues in Root Finding How to Guarantee Robustness? • Rational / Exact Arithmetic • Memory intensive and time consuming • Interval Arithmetic • Inexpensive and robust

41. Robustness Issues in Root Finding Results in Rounded Interval Arithmetic