Subdivision Termination Criteria in Subdivision Multivariate Solvers using Dual Hyper-planes Representations Iddo Hanniel, Gershon Elber Presented by Ofir Markowitz May 2017 1 Abstract - A robust and efficient algorithm for finding solutions (all) of a non-linear multivariate set of equations is highly motivated.

- The simple case of a non-linear univariate equation is already considered a difficult problem, making the multivariate case highly difficult. - A common subdivision-based approach rely on the Bzier or B-spline representation, exploiting their domain clipping and convex hull properties. - After a sub-domain has met some termination criterion, a numerical

improvement method is employed (e.g. Newton-Raphson) or alternatively, the sub-domain is purged due to the criterion determining it contains no solutions. 2 Our Goal Attend the case where the number of constraints/equations is equal to the number of variables - A zero dimensional solution space (discrete solutions).

Present two (main) subdivision termination criteria: 1. The bounding cones criterion: A sufficient condition to determine that a sub-domain has at most one solution. 2. The parallel hyper-planes criterion: The A sufficient condition to determine that a sub-domain has no solutions. 3 Schedule Problem Formal Representation Criterions Control Coefficients Bounding Cones

Parallel Hyper-planes Example Conclusions 4 Problem Formal Representation Given d implicit algebraic constraints: We seek all the solutions simultaneously satisfy them. which

5 Problem Formal Representation We will assume that or B-spline multivariate scalar functions: and are represented as Bzier are the control coefficients. 6 Schedule

Problem Formal Representation Criterions Control Coefficients Bounding Cones Parallel Hyper-planes Example Conclusions 7 Control Coefficients Criterion Firstly, theres a simple termination criterion for determining that a sub-domain has no solutions - check the control

coefficients, i.e., If for some i, the coefficients of are either all positive or all negative, then by the convex hull property of the Bzier/B-spline basis functions, is guaranteed to get only positive or only negative values in the sub-domain, respectively. Thus, one of the constraints is not met and we can drop the subdomain. 8 Schedule Problem Formal Representation Criterions Control Coefficients

Bounding Cones Parallel Hyper-planes Example Conclusions 9 Bounding Cones Criterion 10 Motivation We would like a criterion which determines that a sub-domain has at most one solution.

Lets examine our set of constraints again: Each constraint represent a d-dimensional hyper-surface. So essentially we are looking for the intersection of d of these hypersurfaces. A general hyper-surface is difficult to compute, let alone intersect 11 Cones Definition Let vector

Furthermore, we define the set of vectors on or in the cone: denote the cone with the axis in the direction of a unit-length and an opening angle : Similarly, we define the set of vectors on or outside the cone: 12 Normal & Tangent cones - Definition

For each implicit hyper-surface we denote as its normal cone, i.e., a cone containing the set of all possible normal vectors, and their scalar multiples. Similarly, we denote the complementary cone

(or tangent cone), as the set of all the vectors which are perpendicular to some vector in the normal cone: 13 Normal & Tangent cones - Computation In contrast to an implicit hyper-surface compute. Consider the gradient: its normal cone is easy to Where the terms

denote the coefficients of the partial derivative of with respect to elevated to a common subspace. This is actually a Bzier/Bspline surface with a set of control points in 14 Normal & Tangent cones - Computation By the convex hull property, the set of points in is contained in the control grid generated by the control points: Thus, the computation of a normal bounding cone can be done by: 1. Taking the average of the control vectors (points) in normalized form to

get the cones axis vector . 2. Taking the maximum angle between . Once we computed the normal cone tangent cone is simply done by this: and the control vectors to be , computing the 15

Normal & Tangent cones - Computation The bounding cones obtained by the procedure described here is generally not optimal. There is an algorithm for getting an optimal bounding cone proposed by Gill Barequet and Gershon Elber in The article: Optimal bounding cones of vectors in three and higher dimensions, Information Processing Letters 93 (2005) 8389 We will not discuss it here. 16 Bounding the Surface using the Tangent Cone Now we show that the tangent cone, translated to a point on the surface

bounds the entire surface. 17 Theorem 1 The corollary derived from the claim gives us the following termination criterion: 18 Theorem 1 - Proof 19

Identifying Intersection of Bounding Cones The conditioned derived by Theorem 1 theoretically enables us to determine that a sub-domain has at most one solution. However, we need to compute the intersection of d cones in , which is a difficult problem.

For our purposes, it is sufficient to determine whetherthe tangent cones intersection is trivial (i.e. contains only the origin). 20 New Representation of the Tangent Cones Instead of the (vector, angle) pair representation of the cones, we consider the two hyper-planes intersecting the unit hyper-sphere and the tangent cone.

21 New Representation of the Cones Firstly, well demonstrate how we can determine whether the tangent cones trivially intersect in a 2-dimensional example: 22 New Representation of the Cones We can see that the cones trivially intersect if and only if the parallelogram bounded by the four lines (hyper-planes) is confined to the unit circle. In (b) the cones trivially intersect, and in (c) the cones has a non-trivial intersection. We would like to expand this result to the general case in

and prove it. 23 Cone Representation - Formal Definition Let be a given complementary tangent cone and let denote the unit hyper-sphere in , i.e. A vector on the intersection of i.e.,

and satisfies 24 Cone Representation - Formal Definition Thus, the intersection of the unit hyper-sphere and the tangent cone is delimited between two parallel hyper-planes ] and

Moreover, the intersection of the unit hyper-sphere and the tangent cone is equal to the intersection of the unit hyper-sphere and the space delimited between the two hyper-planes. We will use this fact to our advantage. 25 Cone Representation - Formal Definition Actually, theres a one to one correspondence between dual hyperplanes intersecting the unit hyper-sphere and a tangent cone. Given a dual hyper-planes vector is normalized and tangent cone is

where the coefficient , the corresponding Thus, the dual hyper-planes representation is unique. 26 Intersecting the Tangent Cone and the Unit Hypersphere We want to properly define the space between the two delimiting hyper-planes. Let be the half-space bounded by and containing the

origin, and similarly, let be the half-space bounded by and containing the origin. Note: With ~ : Its a plane. Without ~ : Its a half-space. 27 Intersecting the Tangent Cone and the Unit Hypersphere

Now we can get the intersection between the cone and the unit hypersphere in this fashion: This is a d-dimensional strip on the unit hyper-sphere. 28 Intersecting the Strips The intersection of the d strips yields the following: Hence, the intersection of the strips is actually the intersection of a convex polytope in the unit Hyper-sphere. Well show that all we need to know is

whether theres a vertex of the polytope and 29 Lemma 1 Reminder: 30 Combining Lemma 1 & Theorem 1

Conclusion of Lemma 1 and Theorem 1: There is at most one solution in the sub-domain if and only if the convex polytope is fully confined in the unit hyper-sphere. 31 Combining Lemma 1 & Theorem 1 Consider the degenerate case where the tangent cones axes are linearly dependent. The polytope will be unbounded and thus its intersection with the sphere will

not be empty. If the axes are linearly independent we need to determine whether theres a vertex of the polytope outside the unit hyper-ball. 32 Polytope - Finding the Vertices Each vertex of a d-dimensional polytope is created by simultaneous intersection of d hyper-planes.

The number of possible vertices of a 2d hyper-planes in is . 33 Polytope - Finding the Vertices Luckily, our 2d hyper-planes are

arranged in parallel pairs, so only vertices are possible (It is actually a hyper-parallelepiped). Furthermore, because of the symmetry of our polytope around the origin, if is a vertex then so is and they are the same distance from the origin. And so it suffices to check only half of the vertices, i.e.

34 Linear System of Equations Computing the intersection of d hyperplanes amounts to solving a set of linear equations. Note that these equations are linearly independent because this is not the degenerate case.

Each combination of plus or minus signs generates a vertex, and by the previous discussion, we can solve only half of the combinations, i.e. . 35 Linear System of Equations The coefficient matrix remains

the same for each solution vector, so we can perform QR/LU factorization (or alternatively, matrix inversion) only once. Well assume the factorization/inversion takes operations where 2 < c < 3, and that backward/forward substitution take operations.

Then, computing all the vertices will take operations. 36 Linear System of Equations - Improvement The algorithm for solving the sets of equations can be improved.

Then the previous system of equations can be rewritten as Hence, after inverting V, the solution is: Denote

Thus, the all the for each i. solutions are one of the combinations: 37 Linear System of Equations - Improvement We can encode each solution by d-1 digit binary number where the ones represent a plus sign and the zeros represent a minus sign.

Then, using grey code, we can assure that each iteration, the solution vector only changes by one component. Thus, we can simply add (or subtract) which has changed. Adding (subtracting) a d-tuple takes operations to

, where i is the component , thus we reduce our number of 38 Summarization We cam now summarize the full algorithm for identifying whether a subdomain contains at most one solution: 1. Compute the normal cones,

2. For each cone, extract the pair of delimiting hyper-planes. 3. For each of the a) b) 4. for each hyper-surface.

sets of linear equations described: Solve the set of equations. Let be the ddimensional solution vector. If return False, i.e., the sub-domain is not guaranteed to gave at most a single solution. Return True. The sub-domain is guaranteed to have at most a single solution. 39 Summarization The loop computing the polytope vertices can terminate early if one vertex is out of the unit hyper-sphere.

The most time consuming operation is computing the normal cones, so in order to optimize this process we can compute the gradients once and treat the pair as a (d+1)-dimensional hyper-surface (elevating the partial derivatives to a common function space as ), and subdividing their domain together. 40 Schedule Problem Formal Representation Criterions Control Coefficients

Bounding Cones Parallel Hyper-planes Examples Conclusions 41 Parallel Hyper-planes Criterion 42 Motivation We introduced a criterion for for determining a sub-domain has at

most one root. After this criterion is met, we can apply a numerical improvement step such as Newton-Raphson. But what if the sub-domain has no roots? Then the numerical step may perform a significant number of iterations. 43 Motivation We want to try bounding the solution to an region, and check whether this region is entirely outside the sub-domain. A sub-domain outside the region will be purged. 44

Computing the Graphs Recall the scalar Bzier/B-spline representation of . Their graph is a (d+1)-dimensional hyper-surfaces. Denote their graph as . We can find the control points generating the

hyper-surface in Bzier/B-spline form, using the nodal points (Greville abssica) of the representation of as Bzier/B-spline, i.e. promoting to its graph. 45 Two Bounding Hyper-planes We want to bound the graph

hyper-planes. by two More precisely, we want to bound the region of the sub-domain where . 46 Two Bounding Hyper-planes These hyper-planes are constructed as follows:

47 Two Bounding Hyper-planes These hyper-planes are constructed as follows: 1. Compute the unit normal the center of the sub-domain. of 48 Two Bounding Hyper-planes

These hyper-planes are constructed as follows: 1. Compute the unit normal of the center of the sub-domain. 2. Project the control points onto the normal . 49 Two Bounding Hyper-planes These hyper-planes are constructed as follows: 1. Compute the unit normal

of the center of the sub-domain. 2. Project the control points onto the normal 3. Intersect two (d+1)-dimensional hyperplanes, perpendicular to the normal , through the extreme projected points. 50 Two Bounding Hyper-planes These hyper-planes are constructed as follows: 1. Compute the unit normal of

the center of the sub-domain. 2. Project the control points onto the normal 3. Intersect two (d+1)-dimensional hyperplanes, perpendicular to the normal , through the extreme projected points. 4. Intersect the (d+1)-dimensional hyper-planes with , constructing two ddimensional parallel hyper-planes, denoted 51 Hyper-planes - Formal Construction Denote the point that is maximal (resp. minimal) of

the projected points. Then the parallel hyper-planes are given by 52 Hyper-planes - Formal Construction Finally, the hyper-planes are given by forcing . Then, we have: where 53

The Region Bounding the Solution We would like to define properly the region delimited between and where the solutions is bounded. Denote and as the half-spaces bounded generated by and respectively, such that the area where is in their positive sides. Then, the region containing the

given by is 54 The Region Bounding the Solution We want to intersect the bounding regions of all the constraints so the solutions will be bounded. This intersection is the following: If its not empty, its a d-dimensional polytope (more precisely a hyperparallelepiped). We want to check if its entirely outside the

55 The Region Bounding the Solution There are a number of well-known numerical methods to check this condition. But this problem is similar to the one weve faced before, i.e. checking if one of a hyperparallelepipeds vertices are outside some region, this time its our sub-domain instead of the unit hyper-sphere. Computing the vertices of this polytope will be done in exactly the same fashion. 56

Schedule Problem Formal Representation Criterions Control Coefficients Bounding Cones Parallel Hyper-planes Example Conclusions 57 Example The algorithms was performed on two bivariate Bzier as constraints. The constraints are represented as bicubic Bziers:

The subdivision tolerance is . 58 Example This is the resulting subdivision trees. The black dots are the solutions and the gray dots marks the areas where the subdivision process was terminated due to the bounding cone test. Without the two criteria With the bounding cone criterion

With both criteria 59 Example We zoomed in the black bordered squares where the two close roots resides. We can see a significant improvement when using the bounding cone criterion. Also that the no solution criterion helps in some difficult cases where the bounding cone criterion fails to do so. Without the two criteria With the bounding cone criterion

With both criteria 60 Schedule Problem Formal Representation Criterions Control Coefficients Bounding Cones Parallel Hyper-planes Example Conclusions 61

Conclusions We have presented an efficient and simple algorithm for identifying no or single solution sub-domains in solving a set of d implicit constraints in d unknowns. The algorithms was based on a dual representation of the tangent cones as hyper-planes over the unit hyper-sphere, and on pairs of bounding hyperplanes .

The single solution test guarantees that all the separable roots in the solution set will be isolated correctly. The algorithm offers ways to improve the performance of subdivision based 62