Survey%20and%20Recent%20Results:%20Robust%20Geometric%20Computation

1
Survey and Recent Results Robust Geometric
Computation

• Chee Yap
• Department of Computer Science
• Courant Institute
• New York University

2
OVERVIEW

• Part I NonRobustness Survey
• Part II Exact Geometric Computation
• Core Library
• Constructive Root Bounds
• Part III New Directions
• Moores Law and NonRobustness
• Certification Paradigm
• Conclusion

3
Numerical Nonrobustness Phenomenon
4
Part I OVERVIEW

• The Phenomenon
• What is Geometric?
• Taxonomy of Approaches
• EGC and relatives

5
Numerical Non-robustness

• Non-robustness phenomenon
• crash, inconsistent state, intermittent
• Round-off errors
• benign vs. catastrophic
• quantitative vs. qualitative

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

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 couldnt do it, so we dont have
  to bother

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

9
What is Special about Geometry?
10
Geometry is Harder

• Geometry CombinatoricsNumerics
• E.g. Voronoi Diagram

11
Example Convex Hulls
2
1
2
1
2
1
3
3
3
7
7
8
7
8
4
4
4
9
9
6
6
6
5
5
5
Input
Convex Hull
Output
12
Consistency

• Geometric Object
• Consistency Relation (P)
• E.g. D is convex hull or Voronoi diagram
• Qualitative error ? inconsistency

D (G, L,P) Ggraph, Llabeling of G
13
Examples/Nonexamples

• Consistency is critical
• matrix multiplication
• shortest paths in graphs (e.g. Djikstras
  algorithm)
• sorting and geometric sorting
• Euclidean shortest paths

14
Taxonomy of Approaches
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

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-Stolfi89)
• finite resolution geometry (e.g., grids)
• topology oriented approach Sugihara-Iri88

17
Example

• Grid Geometry Greene-Yao86
• Finite Resolution Geometries

18
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

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?

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 Yu92, Karasick-Lieber-Nac
  kman89

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

, , ?, ?, ?
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)

23
Consistency Approach

• Goal ensure that no decisions are contradictory
• Parsimonious Algorithms Fortune89 only make
  tests that are independent of previous results
• NP-hard but in PSPACE

24
Consistency is Hard

• Hoffmann, Hopcroft, Karasick88
• 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)

???
25
Exact Geometric Computation
26
Part II OVERVIEW

• Exact Geometric Computing (EGC)
• The Core Library
• Constructive Root Bounds

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!

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-Yap94
• 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!

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
  like x1?? x2).
• An expression E? ???? is a DAG
• E ?????????????? with ?? and ?? shared
• Value function, Val ???? ?? R where Val(E)
  may be undefined

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

31
Complexity of CSP

• Hierarchy
• ?0 ? ??, ????? ? ??????????????
• ?1 ? ?0 ? ? ? ?
• ?2 ? ?1 ? ?? ?
• ?3 ? ?2 ? ?Root(P,i) P(x)?Zx ?
• ?4 ? ?0 ? ?Sin(x), ? ?
• Theorem CSP(?3) is decidable.
• Theorem CSP(?1 ) is alternating polynomial time.
• Is CSP(?4) is decidable?

32
Root Bound

• A root bound for an expression E is a
  value such that
• E.g., the Cauchys bound says that
  because is a root of the
  polynomial x4 ???0 x2 ? 1.
• Root bit-bound is defined as ? log(b)

b gt 0
E ¹ 0 Þ E ? b
??????????????
33
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
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.

35
Core Library
36
EGC Libraries

• GOAL any programmer may routinely construct
  robust programs
• Current Libraries
• Real/Expr Yap-Dube94
• LEDA real Burnikel et al99
• Core Library Karamcheti et al99

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

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)

39
Delivery System

• No change of programmer behavior
• At the flip of a switch!
• Benefits code logic verification, fast debugging

define Level N // N1,2,3,4 include
CORE.h /
any C/C Program Here
/
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

41
What is in Level III?

• Fundamental gap between Levels II and III
• Need for iteration consider
• a b c
• Precision sensitive evaluation

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.

43
Precision-Driven Eval of Expressions

• Exprs 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

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.

45
Example Theorem Proving Application

• Kahans example (4/26/00)
• To show that you need theorem proving, or why
  significance arithmetic is doomed
• F(z) if (z0) return 1 else (exp(z)-1)/z
• Q(y) y-sqrt(y2 1) - 1/(y-sqrt(y21))
• G(x) F(Q(x)2).
• Compute G(x) for x15 to 9999.
• Theorem proving with Core Library
  Tulone-Yap-Li00
• Generalized Schwartz Lemma for radical expressions

46
Constructive Root Bounds
47
Problem of Constructive Root Bounds

• Classical root bounds (e.g. Cauchys) 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

48
Illustration

• ? is root of A(X) ?i0 ai Xi of degree n
• Height of A(X) is A ?
• Degree-height Yap-Dube 95
• Cauchys bound
• maintain bound on heights
• but this requires the maintenance of bounds on
  degrees
• Degree-length Yap-Dube 95
• Landaus bound

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

50
BFMS Bound

• Burnikel-Fleischer-Mehlhorn-Schirra99
• For radical expressions
• Tight for division-free expressions
• For those with divisions
• where E is transformed to U(E)/L(E)
• Improvement 01

51
New Constructive Bound

• Applies to general - expressions.
• The bounding function is
• an upper bound on the absolute value
  of conjugates of
• an upper bound on the leading
  coefficient of
• an upper bound on the degree of

52
Inductive Rules
53
Comparative Study

• Major open question is there a root bit-bound
  that depends linearly on degree?
• No single constructive root bound is always
  better than the others.
• BFMS, degree-measure and our bound
• Compare their behavior on interesting classes of
  expressions
• sum of square roots,
• continued fraction, etc.
• In Core Library, we compute all three bounds and
  choose the best one.

54
Experimental Results

• A critical predicate in Fortunes sweepline
  algorithm for Voronoi diagrams is
  .
• Input coordinates are L-bit long, as, bs and
  ds are3L-, 6L- and 2L-bit integers,
  respectively.
• The BFMS bound (79 L 30) bits
• The degree-measure bound (64 L 12) bits
• Our new bound (19 L 9) bits
• Best possible Sellen-Yap (15 L O(1)) bits
  
  

55
Timing on Synthetic Input

• Timings for Fortunes predicate on degenerate
  inputs (in seconds)
• Tested on a Sun UltraSparc (440MHz, 512 MB)

56
Timing for Degenerate input

• Timings on Fortunes algorithm on degenerate
  inputs on a uniform (32 x 32) grid with L-bit
  coordinates
• (in seconds)

57
Moores Law and Non-robustness Trends
58
Part III OVERVIEW

• New Directions
• Robustness as Resource
• (or, how to exploit Moores Law)
• Certification Paradigm
• (or, how to use incorrect algorithms)

59
Moores Law and Robustness

• Computers are becoming faster (Moores Law, 1965)
• Software are becoming less robust
• crashes more often
• must solves larger/harder problems (e.g. Meshes,
  CAD)
• expectation is increasing
• Inverse correlation
• is this inevitable?

60
Reversing the Trend

• Robustness all-or-nothing ?
• Computational Paradigm
• robustness as computational resource
• exchange some speed for robustness
• Goal a positive correlation between Moores Law
  and robustness

61
Robustness Trade-off Curves

• Each program P, for given input, defines a family
  of speed-robustness trade-off curves, one curve
  for each CPU speed

• What is 100 robustness?

62
Computing on the Curve

• P only needs to be recompiled to operate at any
  chosen point on a curve.
• With each new hardware upgrade, we can
  automatically recompile program to achieve the
  sweetspot on the new curve

63
Architecture

• Automatically generated Makefiles
• a registry of programs that wants to join
• sample inputs (various sizes) for each program
• robustness targets/parameters
• optimization function (e.g., min acceptable speed)

64
Hardware Change/Upgrade

• For given hardware configuration, for each
  program
• compile different robust versions of the program
• run each on the test inputs to determine
  trade-off curve
• optimize the parameters
• Layer above the usual compiler
• Application libraries, CPU upgrades
• Can be automatic and transparent (like the
  effects of Moores law)

65
Certification Paradigm
66
Another Paradigm Shift

• Motivation programs based on machine arithmetic
  are fast, but not always correct.
• Floating point filter phenomenon
• Need new framework for working with useful but
  incorrect algorithms

67
Verification vs. Checking

• Program verification movement (1970s)
• Another reason why it failed
• Algorithm always correct
• H-algorithm always fast, often correct
• Checking paradigm Blum,Kannan
• use H-algorithms
• only check specific input/output pairs, not
  programs.

68
Checking vs. Certifying

• Checker given an input-output pair, always
  accept or reject correctly
• Certifier (or Filter) given an input-output
  pair, either accept or unsure
• Why certify?
• Certifiers are easier to design
• Nontrivial checkers may not be known (e.g. sign
  of determinant)

69
Filtered Algorithms

• How to use Certifiers? Ingredients
• Algorithm A
• H-algorithm H
• Filter or Certifier F
• Basic Architecture
• run H, filter with F, run A if necessary

70
Some Floating Point Filters

• FvW Filter, BFS Filter
• Static, Semi-static and dynamic filters
• Lemma
• given expression E with k operations and floating
  point inputs, can compute an upper bound on E
  with 3k flops
• Static case FwW about 10 flops per arithmetic
  operations

71
Model for Filtering

• Case 2D Delaunay Triangulation
• Base line
• Top line
• sigma top line/base line 60
• phi filtering performance 20
• efficacy sigma/phi 3
• Synthetic vs. Realistic algorithms
• beta factor
• filterable fraction of work at top line
• estimating beta Mehlhorn et al

72
Extensions

• Different levels of granularity
• Checking Geometric structures Mehlhorn, et al
• Determinant filter Pan
• Filter compiler Schirra
• Cascaded bank of filter Funke et al
• skip A if necessary

73
Conclusion
74
REVIEW

• Part I NonRobustness Survey
• Part II Exact Geometric Computation
• Core Library
• Constructive Root Bounds
• Part III New Directions
• NonRobustness as Resource (or, how to exploit
  Moores Law)
• Certification Paradigm (or, how to use incorrect
  algorithms)

75
Summary

• Non-robustness has major economic/productivity
  impact
• Scientifically sound and systematic solution
  based on EGC is feasible
• Main open problem
• Is EGC for possible for non-algebraic problems?
  (CSP)
• Program filtering goes beyond program checking
• Robustness as resource paradigm can be exploited
  (with Moores Law)

76
Download Software

• Core Library Version 1.4
• small and easy-to-use
• Project Homepage
• http//cs.nyu.edu/exact/

77
Last Slide
78
Extras
79
Algebraic Degree Bound

• The expression where
  is the radical or root node in the expression.
• The degree bound , where
  is either the index of radical nodes or the
  degree in polynomial root nodes.

80
Leading Coefficients and Conjugates Bound

• Basic tool resultant calculus
• Case given two algebraic numbers and
  with minimal polynomials and
  and a defining polynomial for is
  
• Represented in the Sylvester matrix form
• Leading coefficient
• Constant term
• Degree

81
(cont.)

• Due to the admission of divisions, we also need
  to bound
• Tail coefficients
• Lower bounds on conjugates
• Measures are involved in bounding tail
  coefficients.

82
End of Extra Slides