Un solveur de syst` emes non lin eaires bas e sur le polytope de - - PowerPoint PPT Presentation

▶
un solveur de syst emes non lin eaires bas e sur le
SMART_READER_LITE
LIVE PREVIEW

Un solveur de syst` emes non lin eaires bas e sur le polytope de - - PowerPoint PPT Presentation

Mar 04, 2024 126 likes •517 views

Nonlinear Systems Solver based on a Bernstein Polytope Chr. F unfzig Un solveur de syst` emes non lin eaires bas e sur le polytope de Bernstein Christoph F unfzig Dominique Michelucci Sebti Foufou LE2I, Universit e de

slide-1
SLIDE 1

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Un solveur de syst` emes non lin´ eaires bas´ e sur le polytope de Bernstein

Christoph F¨ unfzig Dominique Michelucci Sebti Foufou LE2I, Universit´ e de Bourgogne, Dijon, France October 2009

1 / 38

slide-2
SLIDE 2

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Overview

Motivation Multivariate Polynomials in the Tensorial Bernstein Basis (TBB) Bernstein Polytope bounding Monomials of Canonical Basis Subdivision Solver using the Bernstein Polytope Handling of LP Solution Errors Examples Conclusion and Future Work

2 / 38

slide-3
SLIDE 3

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Cone Sections (in 2D)

2 variables 2 equations

3 / 38

slide-4
SLIDE 4

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Gough-Stewart Platform

Forward kinematics: Cartesian coords: 3 · 3 variables (3 unknown points, 3 fixed), 9 equations Cayley-Menger determinants: 2 unknown distances d24 and d35, 2 degree-4 equations

4 / 38

slide-5
SLIDE 5

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Tensorial Bernstein-based Solvers

IPP [Sherbrooke etal 1993], Synaps [Mourrain,Pavone 2005]: Similar to the tensorial product for the canonical basis (1, x1, x2

1, . . . xd1 1 ) × (1, x2, x2 2, . . . , xd2 2 ) × . . .

use the tensorial product for Bernstein polynomials (B(d1) (x1), . . . , B(d1)

d1 (x1)) × (B(d2)

(x2), . . . , B(d2)

d2 (x2)) × . . . ◮ Convex hull property of the Bernstein polynomials

extends to multivariate polynomials!

5 / 38

slide-6
SLIDE 6

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Tensorial Bernstein-based Solvers

IPP [Sherbrooke etal 1993], Synaps [Mourrain,Pavone 2005]: Similar to the tensorial product for the canonical basis (1, x1, x2

1, . . . xd1 1 ) × (1, x2, x2 2, . . . , xd2 2 ) × . . .

use the tensorial product for Bernstein polynomials (B(d1) (x1), . . . , B(d1)

d1 (x1)) × (B(d2)

(x2), . . . , B(d2)

d2 (x2)) × . . . ◮ Convex hull property of the Bernstein polynomials

extends to multivariate polynomials!

◮ Current Bernstein-based solvers compute all coefficients

in the TBB! Only the smallest and the largest are relevant for bounding the range or the solution domain!

6 / 38

slide-7
SLIDE 7

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Problem: TBB has exponential cardinality

Polynomials of total degree 2: 1, xi, xixj, x2

i

1 = (B(2)

0 (x1)+. . .+B(2) 2 (x1))×. . .×(B(2) 0 (xn)+. . .+B(2) 2 (xn))

Simplicial Bernstein basis [Reuter et al 2008]: Coverage by and subdivision of simplices! We bound the monomial patches (xi, x2

i ) and (xi, xj, xixj),

i < j by Bernstein polytopes. Polynomial number of faces, still an exponential number of vertices.

7 / 38

slide-8
SLIDE 8

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

The Bernstein polytope for the patch (x, y = x2)

  • ver [0, 1]

1 B0 ≥ 0 B1 ≥ 0 B2 ≥ 0

8 / 38

slide-9
SLIDE 9

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

The Bernstein polytope for the patch (x, y, z = xy) over [0, 1]2

y x z y z x

9 / 38

slide-10
SLIDE 10

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

The Bernstein polytope for the patch (x, y, z = xy) over [0, 1]2

y x z y z x

B(1)

0 (x) × B(1) 0 (y) ≥ 0 → (1 − x)(1 − y) ≥ 0 → 1 − x − y + z ≥ 0

B(1)

0 (x) × B(1) 1 (y) ≥ 0 → (1 − x)y ≥ 0 → y − z ≥ 0

B(1)

1 (x) × B(1) 0 (y) ≥ 0 → x(1 − y) ≥ 0 → x − z ≥ 0

B(1)

1 (x) × B(1) 1 (y) ≥ 0 → xy ≥ 0 → z ≥ 0

10 / 38

slide-11
SLIDE 11

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Our solver

For polynomials of total degree 2, all monomials are : xi, xixj, x2

i . ◮ Associate a variable of an LP to each monomial xi, xixj,

x2

i (called Xi, Xij, Xii).

11 / 38

slide-12
SLIDE 12

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Our solver

For polynomials of total degree 2, all monomials are : xi, xixj, x2

i . ◮ Associate a variable of an LP to each monomial xi, xixj,

x2

i (called Xi, Xij, Xii). ◮ Handle an arbitrary domain i=1...n[ui, vi]:

Scale the Bernstein polytope: Xi = X i−ui

vi−ui with

X i ∈ [ui, vi] and Xi ∈ [0, 1].

12 / 38

slide-13
SLIDE 13

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Our solver

For polynomials of total degree 2, all monomials are : xi, xixj, x2

i . ◮ Associate a variable of an LP to each monomial xi, xixj,

x2

i (called Xi, Xij, Xii). ◮ Handle an arbitrary domain i=1...n[ui, vi]:

Scale the Bernstein polytope: Xi = X i−ui

vi−ui with

X i ∈ [ui, vi] and Xi ∈ [0, 1].

◮ Enclose patches (xi, x2 i ) and (xi, xj, xixj) by their

Bernstein polytopes.

13 / 38

slide-14
SLIDE 14

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Our solver: Box reduction with LP

◮ Each equation and inequality of the system gives a linear

equality or inequality constraint in the LP problem.

14 / 38

slide-15
SLIDE 15

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Our solver: Box reduction with LP

◮ Each equation and inequality of the system gives a linear

equality or inequality constraint in the LP problem.

◮ Reduction: Minimize and maximize every variable Xi.

Some freedom on the order of variables! Reduce in index order once, reduce only largest interval, reduce in equation order, . . .

15 / 38

slide-16
SLIDE 16

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Our solver: Box reduction with LP

◮ Each equation and inequality of the system gives a linear

equality or inequality constraint in the LP problem.

◮ Reduction: Minimize and maximize every variable Xi.

Some freedom on the order of variables! Reduce in index order once, reduce only largest interval, reduce in equation order, . . .

◮ If the LP is not feasible, then the box contains no root.

16 / 38

slide-17
SLIDE 17

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Our solver: Box reduction with LP

◮ Termination: If the domain box is smaller than a

minimum interval width δ in each dimension, then this domain box potentially contains a solution.

17 / 38

slide-18
SLIDE 18

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Our solver: Box reduction with LP

◮ Termination: If the domain box is smaller than a

minimum interval width δ in each dimension, then this domain box potentially contains a solution.

◮ Bisection: If the domain box does not reduce by a

constant factor f = 1

2 then there might be several

  • solutions. Bisect along the longest dimension.

18 / 38

slide-19
SLIDE 19

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Simplex Algorithm in FP Arithmetic

Revised simplex method works through bases defining the LP polytope’s vertices (Vertex coordinates are computed by linear system solves using a LU decomposition of the basis matrix.) Simplex algorithm is affected by FP inaccuracy:

◮ Decision that current basis is optimum!

[Dhiflaoui, Mehlhorn etal SODA 2003]

19 / 38

slide-20
SLIDE 20

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Simplex Algorithm in FP Arithmetic

Revised simplex method works through bases defining the LP polytope’s vertices (Vertex coordinates are computed by linear system solves using a LU decomposition of the basis matrix.) Simplex algorithm is affected by FP inaccuracy:

◮ Decision that current basis is optimum!

[Dhiflaoui, Mehlhorn etal SODA 2003]

◮ Computation of an optimum solution!

20 / 38

slide-21
SLIDE 21

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Handling of LP Solution Errors

Relative error ǫ(x) := |δx|/|x| of the solution vector x by backwards error estimation can be derived using the matrix condition ǫ(x) ≤

κ(A) 1−κ(A)ǫ(A)(ǫ(A) + ǫ(b)) if κ(A)ǫ(A) < 1

where ǫ(A) := |δA|/|A| [Wilkinson 1965] κ(A) := |A||A−1| condition number of A, ǫ(b) := |δb|/|b| rel. error of right-hand side. With the inaccurate borders D(xi) = [[ui, ui]; [vi, vi]]

21 / 38

slide-22
SLIDE 22

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Timings and Examples

22 / 38

slide-23
SLIDE 23

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Circle Packing

Problem

1 2 3

Let G be a completely triangular planar graph, n vertices, 3n − 6 edges. Draw vertices of G as circles. If ij is an edge, then circles i and j are tangent. Circles’ interiors are disjoint. Three points can be fixed arbitrarily for exactly one solution!

◮ If i and j is adjacent:

(xi − xj)2 + (yi − yj)2 − (ri + rj)2 = 0.

23 / 38

slide-24
SLIDE 24

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Circle Packing

Problem

1 2 3

Let G be a completely triangular planar graph, n vertices, 3n − 6 edges. Draw vertices of G as circles. If ij is an edge, then circles i and j are tangent. Circles’ interiors are disjoint. Three points can be fixed arbitrarily for exactly one solution!

◮ If i and j is adjacent:

(xi − xj)2 + (yi − yj)2 − (ri + rj)2 = 0.

◮ If i and j is not adjacent:

(xi − xj)2 + (yi − yj)2 − (ri + rj)2 ≥ 0.

24 / 38

slide-25
SLIDE 25

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Circle Packing

1 2 3 4 5 6 7 8 9 10 11

1 2 1 2

25 / 38

slide-26
SLIDE 26

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Circle Packing: Statistics

Reductions O(n2.2) Effective Reductions O(n1.6) Bisections O(n1.2)

5.2 8.9 12 14.7 20.5 27.4 30.3 43.6 46.9 59.1 65.6 71 78.4 104.5 117.7 172.7 4 6 8 10 12 14 16 18 20 number number of circles Reductions and Bisections (Bernstein inequalities) reductions/10 (reduce all once) effective reductions/10 (reduce all once) bisections (reduce all once) x**2.2 x**1.6 x**1.2

26 / 38

slide-27
SLIDE 27

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Circle Packing: Timings

CPU Timings O(n4.9) with LP code, SoPlex 1.4.1 from ZIB,

  • Berlin. (Intel T7200 2.2 GHz, Windows XP 32Bit)

0.1368 0.7465 4.2108 11.7931 30.6524 75.8022 148.393 232.683 446.952 4 6 8 10 12 14 16 18 20 time [s] number of circles Solution Times Bernstein inequalities (reduce all once) Bernstein inequalities and thick hyperplanes (reduce all once) x**4.9

27 / 38

slide-28
SLIDE 28

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Range Bounds

100 1000 10000 1 2 3 4 5 6 7 8 9 10 Range Width Number of Variables Range Widths with random, zero center-derivative polynomials IA in standard, centered form Bernstein Pn Bernstein P’n Bernstein P’’n TBB TBB Basis TBB Shrinked

28 / 38

slide-29
SLIDE 29

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Range Bounds

10 100 1000 1 2 3 4 5 6 7 8 9 10 Range Width Number of Variables Range Widths with random, generic polynomials IA in standard, centered form Bernstein Pn Bernstein P’n Bernstein P’’n TBB TBB Basis TBB Shrinked

29 / 38

slide-30
SLIDE 30

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Conclusion

⊕ It works for a large number of unknowns n and equations m as the polytope is in N = O(n2) space defined by O(n2) + m halfspaces.

30 / 38

slide-31
SLIDE 31

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Conclusion

⊕ It works for a large number of unknowns n and equations m as the polytope is in N = O(n2) space defined by O(n2) + m halfspaces. ⊕ It can handle under-/over-constrained systems (due to LP).

31 / 38

slide-32
SLIDE 32

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Conclusion

⊕ It works for a large number of unknowns n and equations m as the polytope is in N = O(n2) space defined by O(n2) + m halfspaces. ⊕ It can handle under-/over-constrained systems (due to LP). ⊕ It can handle inequalities (due to LP).

32 / 38

slide-33
SLIDE 33

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Conclusion

⊕ It works for a large number of unknowns n and equations m as the polytope is in N = O(n2) space defined by O(n2) + m halfspaces. ⊕ It can handle under-/over-constrained systems (due to LP). ⊕ It can handle inequalities (due to LP). ⊕ It can be extended to higher polynomial degrees. It can be extended to bounded non-polynomial functions.

33 / 38

slide-34
SLIDE 34

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Conclusion

⊕ It works for a large number of unknowns n and equations m as the polytope is in N = O(n2) space defined by O(n2) + m halfspaces. ⊕ It can handle under-/over-constrained systems (due to LP). ⊕ It can handle inequalities (due to LP). ⊕ It can be extended to higher polynomial degrees. It can be extended to bounded non-polynomial functions. ⊖ It is affected by FP inaccuracy in the LP solutions (in comparison to interval methods) but the polytope does not have close vertices!

34 / 38

slide-35
SLIDE 35

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Future Work

◮ The simplex algorithm is not polynomial time in the

worst case, but for this special LP?

35 / 38

slide-36
SLIDE 36

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Future Work

◮ The simplex algorithm is not polynomial time in the

worst case, but for this special LP?

◮ Exploit special form of the Bernstein polytope for a

compact representation of simplex bases (GPU implementation).

36 / 38

slide-37
SLIDE 37

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Future Work

◮ The simplex algorithm is not polynomial time in the

worst case, but for this special LP?

◮ Exploit special form of the Bernstein polytope for a

compact representation of simplex bases (GPU implementation).

◮ Find good start basis for this special LP: choose linear

independent system rows and fill-up with Bernstein inequalities (e.g., for system rows which contain only x2

i

but not xi)!

37 / 38

slide-38
SLIDE 38

Nonlinear Systems Solver based on a Bernstein Polytope

  • Chr. F¨

unfzig

Thank you!

38 / 38