On Computing Ranges of Polynomials - Some Improvements to the Bernstein Approach

1. On Computing Ranges of Polynomials - Some Improvements to the Bernstein Approach P. S. V. Nataraj Shashwati Ray IIT Bombay Presented at the Second Taylor Model Workshop, Miami, Dec. 2003 Indian Institute of Technology, Bombay 2003

2. SCOPE OF THE PRESENTATION • Introduction • Bernstein forms : • Definition • Properties • Basis conversion • Range calculations • Important theorems • Subdivision • Subdivision direction selection strategies • New propositions • Comparison with Globsol • Examples and conclusion • Future work

3. INTRODUCTION Polynomials are encountered in Engineering and Scientific Applications. Required to find bounds for the range of polynomials over an interval. Our method based on expansion of multivariate polynomial into Bernstein form.

4. ADVANTAGES OF BERNSTEIN FORM • Avoids function evaluations – may be otherwise costly for polynomials of high degrees. • Nice features about range enclosures with Bernstein form: • Information when enclosure is exact. • Sharpness of enclosure can be improved by elevating the Bernstein degree or by subdivision of the domain

5. BERNSTEIN FORM • Any polynomial in power form is represented as • Bernstein form of representation is • Bernstein coefficientsare given by

6. BERNSTEIN FORM • Bernstein functions on [0,1] are defined by • On a general interval the Bernstein functions are defined as • Each set of coefficients ai or bin can be computed from the other (i.e., power form to Bernstein form).

7. BASIS CONVERSION • On unit interval a polynomial’s equivalent power and Bernstein forms • Using matrix multiplication p(x) = XA = BXB • X is the variable (row) matrix. • A is the coefficient (column) matrix • BXis the Bernstein basis (row) matrix • B is the Bernstein coefficient (column) matrix • After certain computations, BXB = XUXB • UX is a lower triangular matrix

8. BASIS CONVERSION • BX = XWXVXUX • WX is an upper triangular matrix • VX is a diagonal matrix • XA = BXB • XA = XWXVXUX B • B=(UX)-1(VX)-1(WX)-1A

9. BASIS CONVERSION • For a bivariate case • X1and X2 are variable matrices • By analogy with univariate case

10. RANGE CALCULATIONIMPORTANT THEOREMS • Theorem 1 : The minimum and maximum Bernstein coefficients give an enclosure of the range of polynomial on the given interval. • Theorem 2 :Vertex condition : If the minimum and maximum Bernstein coefficients of polynomial in Bernstein form occur at the vertices of the Bernstein coefficient array, then the enclosure is the exact range. • Theorem 3 : Bernstein approximations converge to the range and the convergence is at least linear in the order of approximations.

11. SUBDIVISION • Minimum & maximum Bernstein coefficients enclose the range of the polynomial • Vertex condition not satisfied – subdivide in the chosen component direction • Calculate Bernstein coefficient of p(x) on all these subboxes • Check vertex condition on each subbox • Find range of polynomial

12. SWEEP PROCEDURE

13. SUBDIVISION DIRECTION SELECTION STRATEGIESEXISTING RULES • Rule A • Direction chosen is cyclic • Rule B (derivative based) • Find an upper bound for the absolute value of the partial derivative (from its Bernstein form) in the rth direction • Compute the same for each direction • The direction for which this is maximum is selected as the optimal direction for bisection

14. SUBDIVISION DIRECTION SELECTION STRATEGIESEXISTING RULES • Rule C (derivative based) • Find an upper bound for the absolute value of the partial derivative (from its Bernstein form) in the rth direction • Find the product of this upper bound and the degree of rth variable • Compute the same for each direction • The direction for which this is maximum is selected as the optimal direction for bisection

15. SUBDIVISION DIRECTION SELECTION STRATEGIESEXISTING RULES • Rule D (derivative based) • Find an upper bound for the absolute value of the partial derivative (from its Bernstein form) in the rth direction • Find the product of this upper bound and the edge length • Compute the same for each direction • The direction for which this is maximum is selected as the optimal direction for bisection

16. SUBDIVISION DIRECTION SELECTION STRATEGIESEXISTING RULES • Rule E (max width) • Bisection is done along the component direction of maximal width • Rule F (randomized) • Pick a box randomly from the current list • Bisect it in all component directions • Find the maximal width of the polynomial over both these boxes • Choose the component direction where this width is least

17. DISADVANTAGES OF THE EXISTING METHODS OF BERNSTEIN COEFFICIENTS COMPUTATIONS • Conventional method is veryinefficient for higher dimensional polynomials • Even the Matrix multiplication method cannot begeneralized to higher dimensional polynomials

18. DISADVANTAGES OF THE EXISTING METHODS OF BISECTION DIRECTION SELECTION • The existing rules generally need large computation times. • Certain existing rules are grosslyinefficient, especially for higher dimensional polynomials.

19. DISADVANTAGES OF THE EXISTING RANGE COMPUTATION ALGORITHMS • Fortran 95 can not readily cater to more than six dimensional arrays. • Bernstein coefficients generated are stored in multidimensional arrays • For getting sharp enclosures, subdivision usually creates large amount of data • All these are found to heavily slow down computations

20. NEW PROPOSITIONS • We propose new tools for • Bernstein Coefficient Computations • Subdivision Direction Selection • Acceleration devices

21. PROPOSITION 1Computation of Bernstein Coefficients • Matrix Method • This gives

22. PROPOSITION 1Computation of Bernstein coefficients • Extending the above to trivariate case • where transpose means converting 3rd coordinate direction to 2nd , 2nd coordinate direction to 1st , 1st coordinate direction to 3rd • The same logic can be extended to the l-variate case

23. PROPOSED MATRIX METHOD • For a 3-dimensional case instead of considering the polynomial coefficient matrix A as a 3 dimensional array, it can be considered as a matrix with 0 to n1 rows and 0 to (n2+1)(n3+1) –1 columns • After the third transpose and reshape the ‘A’ matrix would come in its original form • The same is shown as figures in the report.

24. Proposed Algorithm Bernstein Matrix • Inputs : • Degree of each variable • Widths of each interval component • Lower end point of each interval component • Polynomial coefficients “a” in matrix form • Output : Bernstein coefficient “a” in matrix form • BEGIN Algorithm • Compute UX, VX and WX for each variable, j =1,2,…, n • Compute their inverses (UX)-1, (VX )-1 and (WX)-1 for j =1,2,…, n • Compute their product Mj for j =1,2,…, n • FOR j =1,2,…, n Set a = Mj a b= transpose(a) a = reshape(b) • END FOR • END Algorithm   

25. Comparison and results (Bernstein coefficients evaluation) • Comparison is done between conventional method and the proposed Matrix method • All problems are taken from Verschelde’s PHC pack (database of polynomial systems) • Both methods generate the same Bernstein coefficients • Proposed method gives results for all the examples • Conventional method takes much longer time; it does not give results for some examples even after a long time.

26. Time in seconds taken by the Matrix and the Conventional methods to compute Bernstein coefficients

27. Some Definitions • Solution box is the box where ‘vertex property’ of the Bernstein coefficients is satisfied. • Solution patch is given by enclosure of the Bernstein coefficients over the solution box. • Current range is the hull of all solution patches.

28. PROPOSITION 2Subdivision Direction Selection • Proposed Rule 1 : Select patch with minimum hull • Find the patches with smallest and largest Bernstein coefficients from the entire list of patches • Find the respective distance of smallest and largest Bernstein coefficient to the infimum and supremum of the current range • Select that patch as the onewhich gives the largest distance, and bisect this selected patch in all component directions

29. EXAMPLEDistance • Let a 3-d polynomial be given as • The polynomial coefficient and the Bernstein coefficient matrices are given as • Let at some stage the Current range be given by [1,20]

30. EXAMPLE • Let the two existing patches in the entire list be • Minimum and maximum value of Bernstein coefficients are 0.5 and 7.70 respectively • The respective distances are (1-0.5) = 0.5 and (20-7.70) = 12.30 • Select the patch which has 7.70 as the Bernstein coefficient

31. Subdivision Direction Selection • Select patch with minimum hull • Hull of the Bernstein coefficients over both the generated subboxes is found out in each direction • The selected bisection direction is the one which gives the minimal hull.

32. PROPOSITION 3Subdivision Direction Selection • Proposed Rule 2 : Randomized box with minimum hull • Choose a box randomly from the current list. • The chosen box is bisected along each and every component direction. • Hull of the Bernstein coefficients over both the generated subboxes is found out in each direction. • The selected bisection direction is the one which gives the minimal hull.

33. RESULTS AND CONCLUSIONSSubdivision Direction Selection • All problems are taken from Verschelde’s PHC pack (database of polynomial systems) • Comparison of the bisection direction selection rules, on the basis of the following performance metrics • Number of solution boxes where vertex property/simplified vertex condition is satisfied • Number of subdivisions • Maximum patch list length • Computation time in seconds • With the proposed rules we are able to solve all the examples • With rules A, B, C, D, E and F we are able to solve 90%, 70%, 70%, 50%, 90% and 90% of the examples, respectively

34. RESULTS AND CONCLUSIONSSubdivision Direction Selection • Comparison is done with all the existing rules (except rule D, which is unable to solve half the problems). • Proposed rules are found to be the most efficient, in terms of every performance metric • Reduction in number of solution boxes compared to Rule A is 82%, Rule B is 83%, Rule C is 83%, Rule E is 81% and Rule F is 84% • Reduction in number of subdivisions compared to Rule A is 84%, Rule B is 83%, Rule C is 83%, Rule E is 81% and Rule F is 84%

35. RESULTS AND CONCLUSIONSSubdivision Direction Selection • Reduction in the maximum list length compared to Rule A is 85%, Rule B is 51%, Rule C is 46%, Rule E is 34% and Rule F is 63% • Reduction in the computational time compared to Rule A is 98%, Rule B is 93%, Rule C is 93%, Rule E is 91% and Rule F is 93% • Both the rules give the same amount of reductions

36. PROPOSITION 4Accelerating Devices for Range Computation • Cut off test • Let the enclosure of Bernstein coefficients over a box be included in the current range. • This box can be deleted from further processing. • Any further subdivision of this box is not going to improve the current range. • Therefore, unnecessary subdivisions are eliminated

37. PROPOSITION 5Accelerating Devices for Range Computation • Simplified Vertex Condition • Let smallest Bernstein coefficient in a patch appear at the vertex ; and the upper bound of this patch be included in the current range • This patch can be taken to be a solution patch • No further subdivision is required, as it would not give any further improvement in the current range • Same is true if the largest Bernstein coefficient in a patch is appearing at the vertex and the lower bound of the patch is included in the current range

38. PROPOSITION 6Accelerating Devices for Range Computation • Monotonicity Test • If the polynomial is monotonic w.r.t any direction on a box • If the box is in the interior • This box can be rejected • If the box is not in the interior • This box is retained • No subdivision of this box in that direction • Avoids unnecessary subdivisions

39. PROPOSED ALGORITHMRange • Inputs : • Maximum degree of each variable of the polynomial • Coefficients of the polynomial in the form of a matrix • Initial domain box • Accuracy tolerances • Outputs: • An enclosure of the range of the specified accuracy • Number of solution boxes • Number of subdivisions • Maximum patch list length used to store the patches

40. PROPOSED ALGORITHMRange • BEGIN Algorithm • Evaluate Bernstein coefficients in matrix form; put them as a patch in a list. • Set flags for simplified vertex condition, cut off test and monotonicity test • Take the first patch from the list and check for ‘vertex’/’simplified vertex’ condition. • If ‘true’, then update current range. • Else, perform ‘cut off test’. • If list is not empty, choose a component direction based on a bisection rule

41. PROPOSED ALGORITHMRange • Take each patch from list, do the monotonicity test. If test is passed, proceed to next step; else discard tested box and choose next patch from the list. • Perform the necessary subdivision in the given direction and evaluate the new Bernstein coefficients • Update the list of patches and remove the tested box • Take the next patch from the list and perform the same operations, till the list is empty • Get the range of the polynomial from the current range • END Algorithm

42. COMPARISON WITH GLOBSOL • We compare our results with the results obtained with Globsol • All problems are taken from Verschelde’s PHC pack, a database of polynomial systems. • Comparison of the computation time required by both the algorithms in terms of values of : • Ratio • Percent reduction

43. COMPARISON WITH GLOBSOL

44. Comparison of computation time by both the algorithms

48. RESULTSComparison with Globsol • In one of the examples Globsol was unable to give any result (arithmetic exception core dumped) indicated by ‘*’ in the table. • Compared to Globsol, the improved Bernstein method gives an average percent reduction in computational time of 47.38%

49. CONCLUSIONS • The new algorithm Range is faster than Globsol in fourteen out of eighteen problems. • Further, on the average, we obtain considerable reductions in computation time with the proposed Algorithm Range.

50. FUTURE WORK • Modify the new algorithm after incorporating certain modifications in the Monotonicity test , thusfurther speeding up the algorithm Range • Integrate the code with COSY package of Berz et al. • Apply the method to Engineering problems, especially control engineering !