Survey and Recent Results: Robust Geometric Computation

1. Survey and Recent Results: Robust Geometric Computation Chee Yap Department of Computer Science Courant Institute New York University Talk @ MIT

2. OVERVIEW • Part I: NonRobustness Survey • Part II: Exact Geometric Computation • Core Library • Constructive Root Bounds • Part III: New Directions • Moore’s Law and NonRobustness • Certification Paradigm • Conclusion Talk @ MIT

3. Numerical Nonrobustness Phenomenon Talk @ MIT

4. Part I: OVERVIEW • The Phenomenon • What is Geometric? • Taxonomy of Approaches • EGC and relatives Talk @ MIT

5. Numerical Non-robustness • Non-robustness phenomenon • crash, inconsistent state, intermittent • Round-off errors • benign vs. catastrophic • quantitative vs. qualitative Talk @ MIT

6. Examples • Intersection of 2 lines • check if intersection point is on line • Mesh Generation • point classification error (dirty meshes) • Trimmed algebraic patches in CAD • bounding curves are approximated leading to topological inconsistencies • Front Tracking Physics Simulation • front surface becomes self-intersecting Talk @ MIT

7. Responses to Non-robustness • “It is a rare event” • “Use more stable algorithms” • “Avoid ill-conditioned inputs” • “Epsilon-tweaking” • “There is no solution” • “Our competitors couldn’t do it, so we don’t have to bother” Talk @ MIT

8. Impact of Non-robustness • Acknowledged, seldom demanded • Economic/productivity Impact • barrier to full automation • scientist/programmer productivity • mission critical computation fail • Patriot missile, Ariadne Rocket • E.g. Mesh generation • a preliminary step for simulations • 1 failure/5 million cells [Aftosmis] • tweak data if failure Talk @ MIT

9. What is Special about Geometry? Talk @ MIT

10. Geometry is Harder • Geometry = Combinatorics+Numerics • E.g. Voronoi Diagram Talk @ MIT

11. Example: Convex Hulls 1 2 1 2 1 2 3 3 3 7 8 7 7 8 4 4 4 9 9 6 6 6 5 5 5 Input Convex Hull Output Talk @ MIT

12. Consistency • Geometric Object • Consistency Relation (P) • E.g. D is convex hull or Voronoi diagram • Qualitative error  inconsistency D = (G, L,P) G=graph, L=labeling of G Talk @ MIT

13. Examples/Nonexamples • Consistency is critical • matrix multiplication • shortest paths in graphs (e.g. Djikstra’s algorithm) • sorting and geometric sorting • Euclidean shortest paths Talk @ MIT

14. Taxonomy of Approaches Talk @ MIT

15. Gold Standard • Must understand the dominant mode of numerical computing • “F.P. Mode” : • machine floating point • fixed precision (single/double/quad) • IEEE 754 Standard • What does the IEEE standard do for nonrobustness? Reduces but not eliminate it. Main contribution is cross-platform predictability. • Historical Note Talk @ MIT

16. Type I Approaches • Basic Philosophy: to make the fast but imprecise (IEEE) arithmetic robust • Taxonomy • arithmetic(FMA, scalar vector, sli, etc) • finite resolution predicates(-tweaking, -predicates [Guibas-Salesin-Stolfi’89]) • finite resolution geometry(e.g., grids) • topology oriented approach[Sugihara-Iri’88] Talk @ MIT

17. Example • Grid Geometry [Greene-Yao’86] • Finite Resolution Geometries Talk @ MIT

18. fat line: polyline: rounded line: aX+bY+c=0; (a<2L, b<2L, c<22L) Example • What is a Finite Resolution Line? • A suitable set of pixels [graphics] • A fat line [generalized intervals] • A polyline [Yao-Greene, Milenkovic, etc] • A rounded line [Sugihara] Talk @ MIT

19. Example • Topology Oriented Approach of Sugihara-Iri:Voronoi diagram of 1 million points • Priority of topological part over numerical part • Identify relevant and maintainable properties: e.g.planarity • Issue: which properties to choose? Talk @ MIT

20. Exact Approach • Idea: avoid all errors • Big number packages (big integers, big rationals, big floats, etc) • Only rational problems are exact • Even this is a problem [Yu’92, Karasick-Lieber-Nackman’89] Talk @ MIT

21. Algebraic Computation • Algebraic number: •     • P(x) = x2 – 2 = 0 • Representation:   P(x), 1, 2) • Exact manipulation: • comparison • arithmetic operations, roots, etc. • Most problems in Computational Geometry textbooks requires only +, –, , ,  Talk @ MIT

22. Type II Approach • Basic Philosophy: to make slow but error-free computation more efficient • Taxonomy • Exact Arithmetic (naïve approach) • Expression packages • Compiler techniques • Consistency approach • Exact Geometric Computation (EGC) Talk @ MIT

23. Consistency Approach • Goal: ensure that no decisions are contradictory • Parsimonious Algorithms [Fortune’89] :only make tests that are independent of previous results • NP-hard but in PSPACE Talk @ MIT

24. Consistency is Hard • [Hoffmann, Hopcroft, Karasick’88] • Geometric Object D = (G, L) • G is realizable if there exists L such that (G, L) is consistent • Algorithm AL:IG(I), L(I)) • AL is geometrically exact if G(I) is the exact structure for input I • AL is consistent if G(I) is realizable • Intersecting 3 convex polygons is hard (geometry theorem proving)  Talk @ MIT

25. Exact Geometric Computation Talk @ MIT

26. Part II: OVERVIEW • Exact Geometric Computing (EGC) • The Core Library • Constructive Root Bounds Talk @ MIT

27. How to Compute Exactly in the Geometric Sense • Algorithm = sequence of steps • construction steps • conditional or branching steps • Branching based on sign   of predicate evaluation • Output combinatorial structure G in D=(G,L) is determined by path • Ensure all branches are correct • this guarantees that G is exact! Talk @ MIT

28. Exact Geometric Computation (EGC) • Exactness in the Geometry, NOT the arithmetic (cf.geometric exactness) • Simple but profound implications • We can now use approximate arithmetic [Dube-Yap’94] • EGC tells us exactly how much precision is needed • No unusual geometries • No need to invent new algorithms -- “standard” algorithms apply • no unusual geometries • General solution (algebraic case) • Not algorithm-specific solutions! Talk @ MIT

29. Constant Expressions •  = set of real algebraic operators. •  = set of expressions over . • E.g., if  , x1 x2 then  is the integer polynomials over x1 x2. • Assume  are constant operators (no variables likex1 x2). • An expression E is a DAG • E = with  and  shared • Value function, Val:   Rwhere Val(E) may be undefined Talk @ MIT

30. Fundamental Problem of EGC • Constant Zero Problem • CZP: given E, is Val(E)=0? • Constant Sign Problem • CSP: given E, compute the sign of Val(E). • CSP is reducible to CZP • Potential exponential gap: • sum of square roots • CZP is polynomial time [Bloemer-Yap] Talk @ MIT

31. Complexity of CSP • Hierarchy: • 0  ,    • 1  0     • 2  1    • 3  2  Root(P,i) : P(x)Z[x]  • 4  0  Sin(x),   • Theorem: CSP(3) is decidable. • Theorem: CSP(1 ) is alternating polynomial time. • Is CSP(4) is decidable? Talk @ MIT

32. Root Bound b > 0 • A root bound for an expression E is a value such that • E.g., the Cauchy’s bound says that because is a root of the polynomial x4 0x2  1. • Root bit-bound is defined as  log(b) E ¹ 0 Þ |E|  b = Talk @ MIT

33.   E E   E   b/2  E How to use root bounds • Let b be a root bound for • Compute a numerical approximation for . • If then sign sign( ) • Else sign (E) = 0. • N.B. root bound is not reached unless sign is really zero! E. ( ) = E Talk @ MIT

34. Nominally Exact Inputs • EGC Inputs are exact and consistent • Why care about exactness if the input is inexact? Because EGC is the easiest method to ensure consistency. Talk @ MIT

35. Core Library Talk @ MIT

36. EGC Libraries • GOAL: any programmer may routinely construct robust programs* • Current Libraries: • Real/Expr [Yap-Dube’94] • LEDA real[Burnikel et al’99] • Core Library [Karamcheti et al’99] Talk @ MIT

37. Core Library • An EGC Library • C++, compact (200 KB) • Focus on “Numerical Core” of EGC • precision sensitive mechanism • automatically incorporates state of art techniques • Key Design Goal: ease of use • “Regular C++ Program” with preamble: #include “CORE.h” • easy to convert existing programs Talk @ MIT

38. Core Accuracy API • Four Levels (I) Machine Accuracy(IEEE standard) (II) Arbitrary Accuracy (e.g. Maple) (III) Guaranteed Accuracy (e.g. Real/Expr) (IV) Mixed Accuracy (for fine control) Talk @ MIT

39. Delivery System • No change of programmer behavior • At the flip of a switch! • Benefits: code logic verification, fast debugging #define Level N // N=1,2,3,4 #include “CORE.h” /* ***************************** * any C/C++ Program Here * ***************************** */ Talk @ MIT

40. What is in CORE levels? • Numerical types: • int, long, float, double • BigInt, BigFloat, Expr • Promotions (+Demotions): • Level II: longBigInt, doubleBigFloat • Level III: long, doubleExpr Talk @ MIT

41. What is in Level III? • Fundamental gap between Levels II and III • Need for iteration: consider a = b + c; • Precision sensitive evaluation Talk @ MIT

42. Relative and Absolute Precision  • Let real X be an approximation of X. • Composite precision bound[r, a] • If r , then get absolute precision a. • If a , then get relative precision r. • Interesting case: [r, a] = [1, ] means we obtain the correct sign of X. Talk @ MIT

43. Precision-Driven Eval of Expressions • Expr’s are DAGs • Each node stores: an approximate BigFloat value; a precision; a root bound • Down-Up Algorithm: • precision p is propagated down • error e propagated upwards • At each node, check e  p • Check passes automatically at leaves • Iterate if check fails; use root bounds to terminate Talk @ MIT

44. Example • Line intersection (2-D): • generate 2500 pairs of lines • compute their intersections • check if intersection lies on one line • 40% failure rate at Level I • In 3-D: • classify pairs of lines as skew, parallel, intersecting, or identical. • At Level I, some pairs are parallel and intersecting, etc. Talk @ MIT

45. Example: Theorem Proving Application • Kahan’s example (4/26/00) • “To show that you need theorem proving, or why significance arithmetic is doomed” • F(z): if (z=0) return 1; else (exp(z)-1)/z; • Q(y): |y-sqrt(y**2 +1)| - 1/(y-sqrt(y**2+1)); • G(x): F(Q(x)**2). • Compute G(x) for x=15 to 9999. • Theorem proving with Core Library [Tulone-Yap-Li’00] • Generalized Schwartz Lemma for radical expressions Talk @ MIT

46. Constructive Root Bounds Talk @ MIT

47. Problem of Constructive Root Bounds • Classical root bounds (e.g. Cauchy’s) are not constructive • Wanted: recursive rules for a family of expressions to maintain parameters p1, p2, etc, for any expression E so that B(p1, p2, ...) is a root bound for E. • We will survey various bounds Talk @ MIT

48.     Illustration •  is root of A(X) = i=0 ai Xi of degree n • Height of A(X) is A  • Degree-height [Yap-Dube 95] • Cauchy’s bound: • maintain bound on heights • but this requires the maintenance of bounds on degrees • Degree-length [Yap-Dube 95] • Landau’s bound: Talk @ MIT

49. Degree-Measure Bounds • [Mignotte 82, Burnikel et al 00] where m( is the measure of . • It turns out that we also need to maintain the degree bound Talk @ MIT

50. BFMS Bound • [Burnikel-Fleischer-Mehlhorn-Schirra’99] • For radical expressions • Tight for division-free expressions • For those with divisions: where E is transformed to U(E)/L(E) • Improvement [‘01] Talk @ MIT