SYMBOLIC PERTURBATION

1. SYMBOLIC PERTURBATION

2. The problem Algorithms in computational geometry make certain assumptions about the input (general status). The treatment of cases violating these assumptions (degeneracies) is tedious, and seldom included in the theoretical discussion. Example: Algorithm constructing convex hull in d dimensions suppose that no d+1 points lie on the same hyperplane.

3. The solution • A perturbation of the input by an infinitesimal factor. • Running the algorithm on the perturb input. • Recover the answer to the original problem from the output of the perturb input.

4. Computational model The input is a set of n vectors in , where and the ith vecot is xi=(xi,1,…,xi,d) for . The basic operations are assumed to be exact. Most algorithms in computational geometry can be modeled by extended algebraic decision tree. The tree is a ternary tree, where each node is labeled by a test function(branch polynomial) and its branches labeled –1,0,+1. Computation is done by evaluating and the branch labeled by the outcome is taken.

5. Computational model – cont. • The set of arithmetic operations computing a branch expression together with the corresponding test if a primitive. • Examples : • Ordering – The comparison of coordinates. The branch polynomial is f = xi,j– xk,j for . • Orientation – The test decides which side of the line a third point lies. f is the determinant of

6. Computational model – cont. • Program – sequence of instructions that implements a specific algorithm. • Execution path – sequence of instructions executed on a particular input. • Complexity measures: • Algebraic model – The total cost is the number of instructions in the longest execution path. • Bit model – The upper bound on the bit complexity of two integers of size O(b) is M(b) = O(b logb loglogb). The total bit complexity is the sum of the costs of every instruction, maximized over all paths.

7. Perturbation The set of degenerate input is represented by and it must be with measure 0. Every input is a point in nd dimensions. The degenerate inputs create a union of surfaces in the nd-dimensional space S with measure 0. For every point x in S, every nonempty open ball around it contains a point y which represents a nondegenerate input, so we can perturb the input that x represents such that all degeneracies disappear. In order not to change any nondegenerate input an open ball small enough is chosen so that it doesn’t intersect any of the above surfaces.

8. Perturbation – cont. A Problem mapping is a mapping between topological spaces. The input space has the standard euclidean topology. The output space Y is (D – finite space(discrete), R – union of real spaces (euclidean)). A degeneracy occurs when the mapping is discontinuous. For any input , a perturbed instance of xis a curve x() rooted at x, i.e the image of a continuous function such that x(0)=x. A perturbation scheme Q defines a perturb instance for every element of x.

9. Perturbation – cont. Definition: Given a problem mapping and a perturbation scheme Q, the perturbed problem mapping is a mapping from X to Y such that assuming that every such limit exists. The goal is that the new problem mapping be defined on a proper superset of the original domain, and that it will be easier than handling the degeneracies. Proposition4: For any perturbation scheme, if mapping Π is continuous at , then .

10. Computing with perturbations • The goal is to obtain from a given program Φ that implements Π, another program that implements . The changes are: • Transforming all arithmetic operations in Φ in order to handle perturbation curves. • Eliminating ε from the output in the postprocessing stage. • Every branching operation f of Φ is transformed to • .

11. Computing with perturbations – cont. Definition: An input instance is degenerate with respect to some program iff it causes some polynomial f at a branch to vanish, where f is not identically zero. Definition: A perturbation scheme is valid with respect to a function f iff, for every input , the limit exists and is nonzero. Perturbation Q is valid with the set of functions iff it is valid with respect to every function in this set. Q is valid with respect to a given program iff it is valid with respect to the set of all branch polynomials in the program.

12. Computing with perturbations – cont. Theorem7: Assume that Q is valid perturbation scheme for a program Φcomputing mapping Πand that is obtained by the transformations described above then (i) computes the perturbed mapping and (ii) for such that Π is continuous, yields Π(x). Statement (ii) holds if some, or all, of the zero branches of Φ are removed. Proof: (i) By validity all limits exist, hence follows the same execution path on x(ε) as Φ would if ε were a small real positive value. Postprocessing is possible so the map is computed by . (ii) is proved by proposition4, and by validity no zero branches are taken in , so they can be pruned away.

13. SOS – Simulation Of Simplicity Orientation: To decide the orientation of a sequence of d+1 points in we use the matrix :

14. SOS-cont. The perturbation is : where for define If 0<ε<1 then (3) Lemma: The set p(ε) is nondegenerate if ε>0 is sufficiently small.

15. SOS-cont. Proof: It is sufficient to show that for no choice of d+1 mutually distinct indices i0,i1,…,id, the determinant of the matrix is equal to zero. We assume and sort the terms of in order of increasing exponents of ε. Specifically, is the first term and the last one.

16. SOS-cont. Each term is of the form bεc. Because ε>0 is arbitrarily small, the absolute value of the first term with nonzero coefficient b is bigger than the sum of all other terms. Such a term exists since (3) guarantees that no two terms have the same exponent of ε, and thus each term cannot be canceled. For example, the coefficient of the last term is not zero, so does not vanish.

17. SOS-cont. Example: d=1. The goal is to calculate the sign of Let k-fold ε-product , ε()=1. Assume i<j Note: the last two terms have no influence on the sign of

18. SOS-cont. If εc= ε((i1,j1),…,(ik,jk)) then we call ε(il,jl) active for The coefficient b can be extracted from the matrix by crossing out all rows and columns that contain an active ε(il,jl). The sign of the coefficient is determined by the number of transpositions needed to sort a certain permutation (j1,j2,…,jd). If it is odd then b=-1 times the determinant. The key observations is that ε((i1,j1),…,(ik,jk)) is the ε-product of a relevant term iff i1<…<ik and j1<…<jk.

19. SOS-cont. The algorithm uses a vector v=[v1,…,vd;vd+1], where each vi corresponds to the ith row of . To encode the ε-product ε((i1,j1),…,(ik,jk)), we set vl=jl. For every i such that the ith row doesn’t contain an active ε(i,j) we define with il the smallest integer in {i1,…,ik,d+1} that is larger than i. Thus, vk implies that ε(k,vk) is active iff vk<vk+1. For example, v=[3,4,4;4] implies that the ε-product of the encoded term is ε(1,3).

20. SOS-cont. In order to generate the terms in the correct order we use the fact that v=[v1,…,vd;vd+1] encodes a more significant term than iff for j, the largest index, such that . this implies that v=[d+1,…,d+1;d+1] encodes the most significant term. Function Next_v(v) returns vector local i,k begin while vi=1 do for down to 1 do return v end

22. Linear Perturbations Linear perturbations are of the following type: xi(ε)=xi+ εbi (2) where and xi=(xi,1,…,xi,d) and bi=(bi,1,…,bi,d) Rn. Let f(x1,…,xn) be any polynomial in n vectors, its initial form is a homogeneous polynomial I(f) equal to the sum of all terms in f of maximum total degree. Theorem8: Let g(x1,…,xn)=I(f) be the initial form of f. For linear perturbation (2) to be valid with respect to f, it suffices that .

23. Linear Perturbations-cont. Proof: Consider f(x(ε)) a univariate polynomail in ε. It is required that f(ε) never vanishes on perturbed input. The highest order coefficient in f(ε) is , which implies that f(ε) has a finite number of roots. It suffices to assume that ε takes real values smaller than the minimum positive root of f(ε). The problem is reduced to finding a single set of input vectors on which the initial form of all the branch polynomials do not vanish. There is also a more generic validity criterion with the perturbation

24. Evaluation of Branch Expressions No computation in evaluating the primitives in perturbed programs need involve the infinitesimal variable. In the case of tests expressed as determinants, the complexity of the evaluation phase dominates the overall complexity. Let MM(t)=O(t2.376) be the algebraic complexity of multiplying two matrices and I the identity matrix. A(ε) is a matrix with polynomial entries in a single variable ε, A(ε) = Arεr + Ar-1εr-1+…+A1ε +A0 where r is the maximum degree in ε of any matrix entry. Ar is always nonsingular.

25. Evaluation-cont. Ar-1A(ε) = Iεr+Ar-1Ar-1εr-1+…+Ar-1A1ε +Ar-1A0 The determinant of the right-hand side equals the characteristic polynomial of Theorem13: Let A(ε) be a matrix of order t, whose entries are linear univariate polynomials in ε; then A(ε)=A1ε+A0, where A0 and A1 are numeric matrices. If A1 is nonsingular, determining the sign of det(A(ε)) can be reduced to computing determinant det(A1) and the characteristic polynomial of matrix –A1-1A0.

26. Evaluation-cont. Proof: By the nonsingularity of A1 we can write A(ε) as A(ε) = (-A1)(-Iε - A1-1A0). If we represent sign with a value in {-1,0,1}, then the sign of det(A(ε)) is the product of det(A1) and det(-Iε-A1-1A0) multiplied by (-1)t. the last determinant is simply the characteristic polynomial of –A1-1A0.

27. Modular Arithmetic Let p0,p1,…,pk-1 be a set of integers that are pairwise relatively prime and let p be the product of the pi. By the Chinese remainder theorem there is one-to-one correspondence between the integers r with and (r0,r1,…,rk-1) with where ri= r mod pi, qi = p/pi, and si=qi-1 mod pi. To evaluate an expression, a set of relatively prime integers is chosen such that the product of the primes it at least twice the absolute value of the integral value of the expression. Then the expression is evaluated modulo each pi. Finally Chinese remaindering is used to reconstruct the value of the expression.

28. Modular Arithmetic – cont. Let k denote the number of primes. The first and third stage each have bit complexity of O(M(k)logk). The middle stage is the actual computation modulo each prime, and its bit complexity is k times the algebraic complexity of this computation. The number of primes is determined dynamically, by choosing random primes and constructing the answer with those primes. At the next iteration a new prime is added and the constructing is done again. If the answer is stable for two or three times then with very high probability the true answer has been calculated.

29. Modular Arithmetic – cont. Corollary14: The algebraic complexity of computing det(A1ε+A0) is O(MM(t)logt), where t is the order of matrices A0 and A1. Let s be the maximum bit size of every entry in A1 and A0. Then the bit complexity of computing the above determinant is , for some arbitrarily small positive constant .

30. Common Primitives • The perturbations dealt with are: • (1) • (5) • Where q is the smallest prime that exceeds n. • Two assumptions in estimating the complexity are made: • Sign determination of determinants is implemented by a determinant calculation. • S=Ω(logn) an upper bound on the bit size of the input.

31. Ordering The ordering primitive decides the order of two quantities expressing the kth coordinate of the i1th and i2th input points. For input perturbed with (1) the primitive decides the sign of The perturbation is valid because all the indices are distinct. The perturbation (5) is valid too for k=1. Theorem: The perturbation defined by (1) is valid with respect to the ordering primitive and doesn’t change the asymptotic running time complexity of this primitive in the algebraic as well as the bit model. The same hold for perturbation (5) for comparisons along the first coordinate.

32. Orientation Given a query point and a hyperplane in Rd spanned by xi1,…,xid, orientation decides in which one of the two half-spaces defined by this hyperplane the query point lies. A degeneracy occurs exactly when lies on the hyperplane. The primitive is formulated as a test of a determinant sign: Proposition: Perturbation (1) is valid with respect to algorithms that branch on determinants of Λδ+1 for , where d is the space dimension. The perturbation increases the asymptotic running-time complexity (algebraic) by O(logd), and the worst-case complexity (bit) by O(d1+α).

33. InSphere InSphere decides, given d+2 points, whether the (d+2)nd point lies in the interior of the higher-dimensional sphere defined by the first d+1 points in Rd. The points are lifted by adding a (d+1)st coordinate equal to the sum of the squares of the d-coordinates defining each point. The lifted images define a hyperplane H in Rd+1. The query point lies within the sphere iff its lifted image lies below H. A degeneracy occurs at the singularities of

34. InSphere-cont. Validity follows by the nonsingularity of matrix Descartes’ Rule Of Sign: The number of sign variations of a polynomial’s nonzero coefficients exceeds the number of positive zeros, multiplicities counted, by an even nonnegative integer. Proposition: Matrix Wd+2 is nonsingular for distinct positive ij, .

35. InSphere-cont. Proof: If Wd+2 is singular, then there is a nonzero vector (q1,…,qd+2) in the kernel of the matrix. Therefore the polynomial has at least d+2 distinct positive zeros, namely i1,…,id+2. The polynomial has at most d+1 sign variations which contradicts Descartes’ rule.

36. InSphere-cont. For the scheme (1), the perturbed determinant expands to Computing each of the two polynomials can be reduced to a characteristic polynomial computation (by theorem13). Theorem: The perturbation of (1) is valid with respect to the InSphere primitive and increases its algebraic complexity by an O(logd) factor. Under the bit model, the worst-case complexity increases by an O(d1+α).

37. Limitations Primitives that decide on the relative position of derived objects may pose a limitation. For example the line-segment intersection problem. The goal is to construct the planar arrangement defined by the segments and their intersections. The hardest primitive is to decide on which half-plane, with respect to a query segment, the intersection of two segments lies.

38. Convex Hull The only tests needed are ordering and orientation, so perturbation (5) is applied to allow arbitrary input. The CHV (convex hull volume) problem is continuous everywhere, so it can be obtained without any postprocessing. The algorithm sorts all points on their first coordinate and then proceeds incrementally by adding each new point to the convex hull of all previous points. Each region between the new point and the existing hull can be partitioned into d-simplices, each defined by the new point and one of the visible facets. The exact volume is computed by summing all simplex volumes.

39. Convex Hull – cont. Perturbation (5) guarantees that the output polytope’s vertex set is a superset of the vertices on the original convex hull. The postprocessing is done by comparing the normals of every to adjacent facets. The normal of a facet can be computed in O(MM(d)) and there are d tests per facet, hence the postprocessing does not affect the asymptotic complexity of the program.