Detecting Degenerate Object Configurations with the Rational Univariate Reduction

1. Detecting Degenerate Object Configurations with the Rational Univariate Reduction John Keyser Department of Computer Science Texas A&M University DIMACS CAD/CAM Workshop

2. Outline • Motivation • The Rational Univariate Reduction • Background • Computation Process • Implementation • Incorporating the RUR • Numerical Perturbation • Conclusion and Future Work DIMACS CAD/CAM Workshop

3. Robustness in Boundary Evaluation • Like many geometric operations, boundary evaluation subject to robustness problems • Exact computation approaches, used in computational geometry, haven’t been applied as readily • Exact computation with curved objects generally requires computer algebra methods. DIMACS CAD/CAM Workshop

4. Goal(big picture) • Support boundary evaluation for curved solids robustly • Numerical error • Degeneracies DIMACS CAD/CAM Workshop

5. Need for Exact Boundary Evaluation • Early work showed numerical problems • Boole system, floating-point based • Curved primitives • Tracing approach, tolerances • Early-mid 90s • Krishnan et al. at UNC DIMACS CAD/CAM Workshop

6. Exact Boundary Evaluation • Solids defined by parametric surfaces • Implicit form also stored • Intersection curves are algebraic plane curves • Usually no parametric representation • Kept in patch domain • Intersections of curves are algebraic points • Solution of systems of polynomial equations DIMACS CAD/CAM Workshop

7. ESOLID • Exact boundary evaluation program • Exact representations of surfaces, curves, points • Key operations: curve-curve intersection, curve topology – provided by MAPC library • MAPC curve-curve intersections • Most time-consuming part of computation • With Culver,Krishnan,Foskey,Manocha DIMACS CAD/CAM Workshop

8. ESOLID (continued) • Intersecting curves • Resultants allow isolating x,y coordinates separately • “Box hits” determine which points are “true” • Points represented as unique root of polynomial system within a box. DIMACS CAD/CAM Workshop

9. ESOLIDSpeedups • Numerous speedup techniques applied • Floating-point filters/Interval arithmetic • Lazy evaluation • Quick rejection • Assumes General Position • Does not handle any actual degeneracies • Fails by crash, infinite loop DIMACS CAD/CAM Workshop

10. ESOLIDPerformance • Applied to real-world cases • Performance, with speedups, within one order of magnitude of Boole DIMACS CAD/CAM Workshop

11. Problem: • Degenerate cases cause failure • Want way to represent points even in degenerate cases • Tangential intersections • Intersections at curve singularities • 3 or more curves meeting at a point • Need way to treat degeneracies DIMACS CAD/CAM Workshop

12. Approach • Define new method for representing points • Allows certain degeneracies to be detected, represented cleanly • Numerical perturbation approach to treat degeneracies • Work done with Koji Ouchi, Maurice Rojas DIMACS CAD/CAM Workshop

13. Outline • Motivation • The Rational Univariate Reduction • Background • Computation Process • Implementation • Incorporating the RUR • Numerical Perturbation • Conclusion and Future Work DIMACS CAD/CAM Workshop

14. Rational Univariate Reduction(RUR) • We will use it as an alternative representation for points • Capable of handling degenerate situations smoothly • Can be used for detecting degeneracies, or as part of a routine to handle them directly DIMACS CAD/CAM Workshop

15. Previous Work on RUR • Fundamental method not new – known for centuries • Prior computation approaches ineffective: • Rely on Groebner bases – too slow • Groebner-free methods are iterative – don’t work well for exact computation • Sparse resultant – Emiris (and others) DIMACS CAD/CAM Workshop

16. Context for Our Implementation • Fits into precision-driven computation model • LEDA and EGC work • Core library • Extends model to handle arbitrary roots of polynomials • Includes complex roots DIMACS CAD/CAM Workshop

17. Outline • Motivation • The Rational Univariate Reduction • Background • Computation Process • Implementation • Incorporating the RUR • Numerical Perturbation • Conclusion and Future Work DIMACS CAD/CAM Workshop

18. RUR operation • Given m polynomials in n variables • Rational coefficients • Determine all roots of system by finding a set of polynomials: h(x): minimal polynomial h1(a), h2(a), etc. coordinate polynomials DIMACS CAD/CAM Workshop

19. RUR operation (continued) • Determine the roots (real and complex) of the minimal polynomial • Evaluate those roots in coordinate equations • Result gives coordinates of every common root of original system • For positive dimensional components, gives one point on that component. DIMACS CAD/CAM Workshop

20. Procedure for Computing RUR Determine Matrix Labels Minimal Polynomial Coordinate Polynomials Compute Perturbed System Input Polynomials Perturb For Each Coordinate DIMACS CAD/CAM Workshop

21. Computing the Matrix Labels • Form and compute the toric perturbation (sparse resultant matrix following Canny and Emiris) • Find Newton polytopes for input equations • Compute Minkowski sum • Perturb randomly • Compute mixed subdivision • Lattice points in interior give column labels • Rows are indexed by monomial multiplied by original equation • Entries are coefficients of input polynomials DIMACS CAD/CAM Workshop

22. Using the Sparse Matrix • Determinant is the sparse resultant • Vanishing is a necessary condition for common root • Determining matrix row/column labels • linear programming • Do not actually put equation coefficients in • Instead, will use perturbed system DIMACS CAD/CAM Workshop

23. Procedure for Computing RUR Determine Matrix Labels Minimal Polynomial Coordinate Polynomials Compute Perturbed System Input Polynomials Perturb For Each Coordinate DIMACS CAD/CAM Workshop

24. Computing the Minimal Polynomial • If input system is f1=f2=…=fn=0 in variables x1,…,xn Compute perturbed system: f0=u0+u1x1+u2x2+…+unxn fi’=fi-sfi* for i=1 to n ui are chosen randomly fi* are polynomials with same support as fi, but with random coefficients DIMACS CAD/CAM Workshop

25. Computing the Minimal Polynomial (continued) • Find sparse resultant of perturbed system • The minimal polynomial is found from this (perturbed) resultant, substiuting a variable T for u0. DIMACS CAD/CAM Workshop

26. Procedure for Computing RUR Determine Matrix Labels Minimal Polynomial Coordinate Polynomials Compute Perturbed System Input Polynomials Perturb For Each Coordinate DIMACS CAD/CAM Workshop

27. Computing Coordinate Polynomials • To compute the coordinate polynomial hi, we substitute ui+/-1 for ui, and compute the determinant again (twice). • We find the square-free parts, and using the two determinants, we again introduce the variable T, with a gcd calculation. • Finally, the coefficient polynomial is found through division. DIMACS CAD/CAM Workshop

28. Computing Polynomials(Summary) • Use linear programming to determine sparse matrix • Evaluate that matrix several times for various entries • Perform some polynomial operations (gcd, division) to determine RUR. DIMACS CAD/CAM Workshop

29. ExtensionsAffine Roots • Can compute affine roots (coordinate is zero) by appending origin to the support of each input polynomial • May lead to spurious additional roots • Will need to test each root found to see whether it should be kept (i.e. solves the actual system). DIMACS CAD/CAM Workshop

30. ExtensionsNon-square Systems • Can deal with non-square systems. • If too few polynomials, just repeat one polynomial • If too many polynomials, form new set of polynomials from generic linear combinations of input polynomials • Will need to test final results, as spurious roots may be found. DIMACS CAD/CAM Workshop

31. Computing with the RUR • Given the RUR, can find roots of minimal polynomial • We need all roots, real and complex • Various techniques – Aberth’s method is one for iteratively converging to roots. • Substitute these into coordinate polynomials to determine (complex) coordinates of roots. DIMACS CAD/CAM Workshop

32. an,x a2,x an,y a2,y an,z a2,z Procedure for Computing RUR Find Roots (Aberth’s Method) Minimal Polynomial … a1 a2 an x-coordinate Polynomial y-coordinate Polynomial z-coordinate Polynomial a1,x a1,y … a1,z DIMACS CAD/CAM Workshop

33. Dealing with Complex Numbers • Do not fit into existing root-bound approaches • Can determine real/imaginary parts separately • Root bounds can be determined for real and imaginary parts independently DIMACS CAD/CAM Workshop

34. Finding Real Roots • Often only care about real roots • To find which roots are real, use root bounds to show that imaginary part is 0. • Remember complex roots come in conjugate pairs • In many cases, no complex roots near 0 – this test is easy. DIMACS CAD/CAM Workshop

35. Representing Roots • The evaluation process can determine bounding intervals (boxes) around each root • Can be used in geometric computations just like previous methods – e.g. for quick rejection tests DIMACS CAD/CAM Workshop

36. Functionality • Each root is found by one point • Degeneracies handled cleanly: • Tangential intersections • Singularities in curves/surfaces • Points lying on curves/surfaces • Coincident points/curves/surfaces DIMACS CAD/CAM Workshop

37. Functionality • At least one point found on each positive dimensional component • Two determinations with different generic values will output two different sets of points • Differences indicate presence of positive dimensional components DIMACS CAD/CAM Workshop

38. Outline • Motivation • The Rational Univariate Reduction • Background • Computation Process • Implementation • Incorporating the RUR • Numerical Perturbation • Conclusion and Future Work DIMACS CAD/CAM Workshop

39. Performance • Implemented in C++, using Gnu MP for multiprecision arithmetic • Exact implementation of sparse resultant • Based off of Emiris’ implementation • Uses exact computation throughout DIMACS CAD/CAM Workshop

40. Determining Polynomials • Use interpolation of polynomials to avoid symbolic computation. • Substitute (random) values and determine • Vandermonde interpolation • Used when evaluating sparse resultant matrix determinant. • Have degree bounds for all polynomials DIMACS CAD/CAM Workshop

41. Timing Results • Quadric curve intersections arising from real-world boundary evaluation cases: • MAPC: .017 -.024 seconds • RUR: .317-1.772 seconds • Approximately 20-100 times slower! • Cases with degenerate intersections, positive dimensional components, higher dimension, all successful DIMACS CAD/CAM Workshop

42. Timing Breakdown • Slows rapidly with higher degrees/dimensions • For lower dimension/degree, the size of the matrix tends to determine running time • Increases quickly with degree/dimension • At higher dimension/degree, coefficient size of coordinate polynomials grows very quickly and tends to dominate time DIMACS CAD/CAM Workshop

43. Checking Root Bounds • Surprisingly, >98% of time was spent in computation of the RUR itself, not in checking root bounds. • Root bounds could conceivably take longer, with repeated construction in precision-driven system. DIMACS CAD/CAM Workshop

44. Optimizations • Other optimizations are notimplemented yet. • Prior experience shows filtering and similar approaches can yield significant speedups • Difficult to filter the sparse resultant matrix calculations • May be able to filter over coefficients of the coefficient equations • Possibly make higher degree/dimension more practical, but unlikely to ever beat MAPC DIMACS CAD/CAM Workshop

45. Outline • Motivation • The Rational Univariate Reduction • Background • Computation Process • Implementation • Incorporating the RUR • Numerical Perturbation • Conclusion and Future Work DIMACS CAD/CAM Workshop

46. Direct Implementation • RUR approach could be used always • Represent all points using RUR • All degeneracies handled smoothly • Checks between points are easy • Evaluating point on curve is easy • However, much slower for cases where there is no degeneracy DIMACS CAD/CAM Workshop

47. Incorporating for Degenerate Situations • Idea: hybrid representation of points • Use standard MAPC approach first • When it does not seem to be working (e.g. repeated iterations), determine RUR • Thereafter RUR can be used instead • Would allow a uniform treatment of all points, with lower overhead than all-RUR • Bounding boxes for both approaches • Might not need RUR itself after first evaluation DIMACS CAD/CAM Workshop

48. Detecting Degeneracies • Several degenerate situations in solid modeling (boundary evaluation) • Some are dealt with more easily by other means: • Overlapping surfaces dealt with through factorization or vanishing of intersection curve • Usually want to detect degeneracies • Apply special case code or perturbation DIMACS CAD/CAM Workshop

49. Degeneracies DIMACS CAD/CAM Workshop

50. Specific Degeneracy Detection • Assume points represented as RUR • Coincident Points: • Check sign (need root bounds) for difference of two points • Four or more surfaces meeting at a point • Appears as a point lying on a curve in a patch domain • Substitute coordinate values of point into curve • Can generate root bounds to guarantee answer DIMACS CAD/CAM Workshop