next up previous contents index
Next: 4.5 Auxiliary variable method Up: 4. Nonlinear Polynomial Solvers Previous: 4.3.3 Subdivision Methods   Contents   Index


4.4 Projected Polyhedron algorithm

In this section we introduce an iterative global root-finding algorithm for an $ n$ -dimensional nonlinear polynomial equation system, which belongs to the class of subdivision methods, called Projected Polyhedron (PP) algorithm developed by Sherbrooke and Patrikalakis [392]. It is easy to visualize and simple in that it only requires two straightforward algorithms in order to implement it: one for subdividing multivariate polynomials in Bernstein form, and one for finding the convex hull of a two-dimensional set of points. This algorithm is an extension and generalization of earlier adaptive subdivision algorithms: for $ n=1$ used in finding the real roots and extrema of a polynomial within an interval by Lane and Riesenfeld [221], and for $ n=2$ used in shape interrogation by Geisow [124] or in intersecting rays with trimmed rational polynomial surface patches by Nishita et al. [286] (a method known as Bézier clipping). The PP algorithm has found many applications in shape interrogation problems (see also Grandine and Klein [133]) as we will see in subsequent sections and its convergence, rate of convergence and complexity properties are developed in [392].

For illustration, we will enumerate the procedures required by the PP algorithm to find roots of a degree $ m$ polynomial equation $ f(x) = c_o + c_1x + c_2x^2 + \cdots + c_mx^m = 0$ over the interval $ a \leq x \leq b$ .

  1. Make an affine parameter transformation $ x = a + t(b - a)$ such that $ 0\leq t \leq 1$ as follows:
    $\displaystyle f(t) = \sum_{i=0}^{m}c_i^Mt^i,\;\;\;\;\;\;0\leq t \leq 1\;.$     (4.9)

    The transition from the interval $ a \leq x \leq b$ to the interval $ 0\leq t \leq 1$ is an affine map, and the polynomials are invariant under affine parameter transformation [92].
  2. Convert the basis from monomial to Bernstein [106]:
    $\displaystyle f(t) = \sum_{i=0}^{m}c_i^BB_{i,m}(t)\;,$     (4.10)

    where
    $\displaystyle c^{B}_{i} = \sum_{j=0}^{i}\frac{(^{i}_{j})}{(^{m}_{j})}c^{M}_{j}\;,$     (4.11)

    and $ B_{i,m}(t)$ is the $ i$ th Bernstein polynomial of degree $ m$ .
  3. Create a graph of function $ f(t)$ using the linear precision property of Bernstein polynomials (see (1.21)). Then the graph will become a Bézier curve
    $\displaystyle {\bf f}(t) = \begin{pmatrix} t \cr f(t) \cr \end{pmatrix} = \sum_{i=0}^{m}\begin{pmatrix}\frac{i}{m} \cr c^B_i \cr \end{pmatrix} B_{i,m}(t)\;,$     (4.12)

    where $ (\frac{i}{m}, c^{B}_{i})^T$ are control points. Now the problem of finding roots of the univariate polynomial has been transformed into a geometric problem of finding the intersection of a Bézier curve with the parameter axis, a transformation already used in Geisow [124] for surface interrogation.
  4. Construct the convex hull of the Bézier curve.
  5. Intersect the convex hull with the parameter axis.
  6. Discard the regions which do not contain roots by applying the de Casteljau subdivision algorithm and find a sub-region of [0,1] which may contain the root(s).
  7. If the sub-region is sufficiently small, we conclude that there is a root inside and return it. But when there are more than one root in the sub-region, the sub-region will not be reduced. In such case we split the region evenly by applying the de Casteljau subdivision algorithm and we go back to 4.

Example 4.4.1. PP algorithm in one dimension.
Find the roots of $ f(x) = -1.1x^2 + 1.4x - 0.2=0$ where $ 0\leq x \leq 2$ . The roots are approximately, 0.164 and 1.108.

  1. Make an affine parameter transformation by plugging $ x = 0 + (2 -0)t = 2t$ into $ f(x) $ yielding $ f(t) = -4.4t^2 + 2.8t -0.2=0$ where $ 0\leq t \leq 1$ .
  2. Convert from monomial to Bernstein basis using (4.11) as below
    $\displaystyle c_i^B = \sum_{j=0}^{i}\frac{(^i_j)}{(^2_j)}c_j^M\;,$      

    where $ c_0^M=-0.2$ , $ c_1^M=2.8$ and $ c_2^M=-4.4$ , thus leading to $ c_0^B=-0.2$ , $ c_1^B=1.2$ and $ c_2^B=-1.8$ .
  3. Create a graph of function $ f(t)$ using linear precision property of the Bernstein polynomial
    $\displaystyle t = \sum_{i=0}^{2}\frac{i}{2}B_{i,2}(t)\;,$      

    yielding a Bézier curve
    $\displaystyle {\bf f}(t) = \begin{pmatrix}t \cr f(t) \cr \end{pmatrix} = \sum_{i=0}^{2}\begin{pmatrix}\frac{i}{2} \cr c^B_i \cr \end{pmatrix} B_{i,2}(t)\;.$      

    with control points of (0, -0.2), (0.5, 1.2) and (1, -1.8).
  4. Construct a convex hull of the Bézier curve, which is a triangle as shown in Fig. 4.3.
  5. The convex hull intersects the t-axis with $ t= 0.0714$ and $ t= 0.7$ .
  6. Discard the regions $ 0\leq t \leq 0.0714$ and $ 0.7\leq t \leq 1$ , which do not contain roots, by applying the de Casteljau algorithm. Now we have a smaller convex hull which contains the roots (see shaded triangular in Fig. 4.3).
  7. If the sub-region is sufficiently small, we conclude that there is a root inside and return it. In this case there are two roots in the convex hull and the sub-region does not reduce much (even after several iteration steps), thus we split the region evenly by applying the de Casteljau subdivision algorithm and go back to 4.
Figure 4.3: de Casteljau algorithm applied to the quadratic Bézier curve
\begin{figure}\centerline{ \psfig{figure=fig/cvx.eps,height=3in} } \end{figure}


Example 4.4.2. PP algorithm in two dimensions.
Let us solve the system of polynomial equations (4.2) and (4.3) over the region $ -1\leq x_1, x_2 \leq 1$ .

  1. Make an affine parameter transformation by substituting $ x_1=2u-1$ and $ x_2=2v-1$ into (4.2) and (4.3) so that $ 0\leq u,v \leq 1$ , then:
    $\displaystyle {f_1}(u,v)=4u^2 - 4u + 4v^2 - 4v + \frac{23}{16} = 0\;,$      
    $\displaystyle {f_2}(u,v)=4u^2 - 8u + 4v^2 - 4v + \frac{19}{4} = 0\;.$      

  2. Convert from monomial to Bernstein basis using
    $\displaystyle c_{ij}^B = \sum_{k=0}^{i}\sum_{l=0}^{j}\frac{(^i_k)(^j_l)}{(^m_k)(^n_l)}c_{kl}^M\;,$      

    where in this case $ m=n=2$ leading to

    \begin{displaymath} \begin{array}{lll} {c_1}_{00}^B = 1.4375,\; & {c_1}_{01}^B =... ...c_1}_{21}^B = -0.5625,\; & {c_1}_{22}^B = 1.4375\;, \end{array}\end{displaymath}


    and

    \begin{displaymath} \begin{array}{lll} {c_2}_{00}^B = 4.75,\; & {c_2}_{01}^B = 2... ... & {c_2}_{21}^B = -1.25,\; & {c_2}_{22}^B = 0.75\;. \end{array}\end{displaymath}

  3. Create graphs of functions $ f_1(u, v)$ and $ f_2(u,v)$ using the linear precision property of Bernstein polynomials. Then the graphs will become two Bézier surfaces as follows:
    $\displaystyle {\bf f}_1(u,v) = \begin{pmatrix}u \cr v \cr {f_1}(u,v) \cr \end{p... ...i}{2} \cr \frac{j}{2} \cr {c_1}^B_{ij} \cr \end{pmatrix}B_{i,2}(u)B_{j,2}(v)\;,$      
    $\displaystyle {\bf f}_2(u,v) = \begin{pmatrix}u \cr v \cr {f_2}(u,v) \cr \end{p... ...}{2} \cr \frac{j}{2} \cr {c_2}^B_{ij} \cr \end{pmatrix} B_{i,2}(u)B_{j,2}(v)\;.$      

    The two Bézier surfaces are shown in Fig. 4.4. Now the root-finding problem of the bivariate polynomial system has been transformed to find the intersections of three surfaces, $ {\bf f}_1(u,v)$ , $ {\bf f}_2(u,v)$ and the $ xy$ -plane. Figure 4.5 shows the intersection between the plane and both Bézier surfaces. We can easily observe that the two intersection curves are the circles in Fig. 4.1 but trimmed in the resulting $ (u,v)$ domain.

    Figure 4.4: Bézier surfaces and their control points
    \begin{figure}\vspace*{15mm} \centerline{ \psfig{file=fig/cir1_bez.eps,height=3.4in} \hspace*{0.2in} \psfig{file=fig/cir2_bez.eps,height=2.87in}}\end{figure}

    Figure 4.5: Bézier surfaces intersecting with $ xy$ -plane
    \begin{figure}\vspace*{0.5in} \centerline{ \psfig{file=fig/cir1_plane.eps,height... ...=2.7in} \psfig{file=fig/cir2_plane.eps,height=2.7in,width=2.7in} }\end{figure}

  4. Project the control points of $ {\bf f}_1(u,v)$ and $ {\bf f}_2(u,v)$ onto $ xz$ and $ yz$ planes. Here $ x=u$ , $ y=v$ . For each $ xz$ and $ yz$ plane, construct 2-D convex hulls. Figure 4.6 (a) shows 2-D convex hulls on the $ xz$ plane, while Fig. 4.6 (b) shows 2-D convex hulls on the $ yz$ plane. The solid line corresponds to convex hull of $ {\bf f}_1(u,v)$ and the dashed line corresponds to that of $ {\bf f}_2(u,v)$ .
    Figure 4.6: Projections of control points
    \begin{figure}\centerline{ \psfig{file=fig/cvh_pro.eps,height=5.0in} }\end{figure}

  5. Intersect each 2-D convex hull on the $ xz$ plane with $ x$ ($ u$ ) axis. The parameter interval $ [0, 1]$ contains solutions of (4.2), while the interval $ [0.34375,1]$ contains solutions of (4.3). The root is contained in the common interval $ [0.34375,1]$ of $ u$ . We will repeat the same procedures for the 2-D convex hulls on the $ yz$ plane to find the common interval $ [0.1875, 0.8125]$ of $ v$ .
  6. Discard the reigion $ [0,0.34375]$ of $ u$ and regions $ [0,0.1875]$ $ [0.8125,1]$ of $ v$ which do not contain the roots by applying the de Casteljau subdivision algorithm to both Bézier surfaces.
  7. If the sub-regions of both parameters are sufficiently small, we conclude that there is a root. In this case, the sub-region of the $ v$ parameter will not decrease much in size because there are two roots, while the interval decreases more in $ u$ parameter, since the two roots have the same $ u$ value. We split the sub-region of $ v$ evenly by applying the de Casteljau subdivision algorithm for both surfaces and we go back to step 4.

Many applications in shape interrogation result in systems of $ n$ nonlinear polynomial equations with $ n$ unknowns, referred to as balanced systems. However, there exist some problems such as tangential intersection or implicit curve/surface rendering consisting of $ n$ nonlinear polynomial equations with $ l$ unknowns, where $ n$ could be larger or smaller than $ l$ . When $ n>l$ the system is called overconstrained and when $ n<l$ it is called underconstrained. Now we will introduce an $ n$ -dimensional Projected-Polyhedron algorithm such that it can effectively handle overconstrained problems [392,179]. This algorithm can also be used in underconstrained problems but in such cases it tends to be slow (especially in the presence of infinite roots); for such cases more specialized algorithms are necessary (e.g. parametric surface intersections where $ n=3$ , and $ l=4$ ). In such cases the PP algorithm is used to find characteristic points efficiently and marching methods are used to trace the intersection curves (see Sect. 5.8.2).

Suppose we are given a set of $ n$ nonlinear polynomials $ f_{1}, f_{2}, \ldots, f_{n}$ , each of which is polynomial in the independent variables $ x_{1}, x_{2}, \ldots, x_{l}$ . Let $ m_{i}^{(k)}$ denote the degree in $ x_{i}$ of polynomial $ f_{k}$ , so that the multi-index $ M^{(k)} = (m_{1}^{(k)}, m_2^{(k)}, \ldots, m_l^{(k)})$ describes all the degree information of $ f_{k}$ . Furthermore, suppose we are given an $ l$ -dimensional rectangular subset of $ {\bf R}^{l}$

$\displaystyle B = [a_{1}, b_{1}] \times [a_{2}, b_{2}] \times \ldots \times [a_{l}, b_{l}]\;.$     (4.13)

A priori knowledge of $ B$ is one of the main features of geometric modeling and shape interrogation problems. We wish to find all points $ {\bf x} = (x_{1}, x_{2}, \ldots, x_{l}) \in B$ such that
$\displaystyle f_{1}({\bf x}) = f_{2}({\bf x}) = \ldots = f_{n}({\bf x}) = 0\;.$     (4.14)

By making the affine parameter transformation [92] $ x_{i} = a_{i} + u_{i} (b_{i} - a_{i})$ for each $ i$ between $ 1$ and $ l$ inclusive, we simplify the problem to one of determining all $ {\bf u} \in [0,1]^{l}$ such that

$\displaystyle f_{1}({\bf u}) = f_{2}({\bf u}) = \ldots = f_{n}({\bf u}) = 0\;.$     (4.15)

Since all of the $ f_{k}$ are polynomial in each of the $ l$ independent variables, a simple change of basis [92] allows us to express them in the multivariate Bernstein basis, which has better numerical stability under perturbation of its coefficients than the power basis [105] as we discussed in Sect. 1.3.3 and in addition permits transformation of an algebraic problem to a geometric problem. In other words, for each $ f_{k}$ there exists an $ l$ -dimensional array of real coefficients $ w_{i_{1}i_{2}\ldots i_{l}}^{(k)}$ such that for each $ k \in \{1, \ldots, n\}$


$\displaystyle f_k ({\bf u}) = \sum_{i_1=0}^{m_1^{(k)}} \sum_{i_2=0}^{m_2^{(k)}}... ...{i_1,m_1^{(k)}} (u_1) B_{i_2,m_2^{(k)}} (u_2) \ldots B_{i_l,m_l^{(k)}} (u_l)\;.$     (4.16)

The notation in (4.16) may be simplified by letting $ I = (i_1, i_2, \ldots, i_l)$ , $ M^{(k)}$ = ($ m_1^{(k)}$ ,$ m_2^{(k)}$ ,$ \ldots$ ,$ m_l^{(k)}$ ) and writing (4.16) in the equivalent form


$\displaystyle f_{k}({\bf u}) = \sum_{I}^{M^{(k)}} w_{I}^{(k)}B_{I,M^{(k)}}({\bf u})\;.$     (4.17)

Representation of algebraic and piecewise algebraic surfaces (i.e., for $ l=3$ ) in terms of tensor products of Bernstein polynomials or B-splines has been studied earlier by Patrikalakis and Kriezis [299]. Equation (4.17) is simply an extension to $ n$ dimensions. Provided that conversion of the problem to the Bernstein basis is exact or sufficiently accurate, the use of the Bernstein basis in conjunction with subdivision is known to be numerically stable [105]. The conversion process itself may be numerically ill-conditioned [106]. Therefore, we recommend that the problem be formulated in the Bernstein basis from the very beginning or that the conversion is carried out in exact arithmetic. If necessary, polynomials may be converted from the multivariate power basis to the multivariate Bernstein basis using the following formula:

$\displaystyle c_{i_1i_2\ldots i_l}^B = \sum_{j_1=0}^{i_1}\sum_{j_2=0}^{i_2}\ldo... ...)}{(^{m_1}_{j_1}) (^{m_2}_{j_2})\ldots (^{m_l}_{j_l})}c_{j_1j_2\ldots j_l}^M\;.$     (4.18)

We now restate the problem as the intersection of the graphs of the $ f_k$ (each of which is a hypersurface in $ {\bf R}^{l+1})$ and the hyperplane $ u_{l+1} = 0$ of $ {\bf R}^{l+1}$ . This idea is designed to impart geometrical significance to the coefficients of the polynomials and to the solution process.

Let us build a graph $ {\bf f}_k$ for each $ f_k$ :

$\displaystyle {\bf f}_k ({\bf u}) = (u_1, u_2, \ldots, u_l, f_k({\bf u}))^T %\nonumber \\ = ({\bf u}, f_k({\bf u}))^T\;.$     (4.19)

Clearly, (4.15) is satisfied by the point $ {\bf u}$ if and only if
$\displaystyle {\bf f}_1 ({\bf u}) = {\bf f}_2 ({\bf u}) = \ldots = {\bf f}_n ({\bf u}) = ({\bf u}, 0)\;.$     (4.20)

Using the linear precision property of the Bernstein basis (1.21), we obtain an equivalent expression for each of the $ u_{j}$ in equation (4.19):

$\displaystyle u_j = \sum_{I}^{M^{(k)}} \frac{i_{j}}{m_j^{(k)}}B_{I,M^{(k)}}({\bf u})\;.$     (4.21)

Substituting (4.21) into (4.19) gives a more useful representation for the $ {\bf f}_k$ :
$\displaystyle {\bf f}_k({\bf u}) = \sum_{I}^{M^{(k)}} {\bf v}_{I}^{(k)}B_{I,M^{(k)}}({\bf u})\;,$     (4.22)

where
$\displaystyle {\bf v}_{I}^{(k)} = \left(\frac{i_{1}}{m_1^{(k)}}, \frac{i_2}{m_2^{(k)}}, \ldots, \frac{i_l}{m_l^{(k)}}, w_I^{(k)}\right)^T\;.$     (4.23)

The $ {\bf v}_I^{(k)}$ are called the control points of $ {\bf f}_k$ . Using the parametric hypersurfaces $ {\bf f}_k$ instead of the real-valued $ f_k$ permits use of the powerful convex-hull property of the multivariate Bernstein basis.

We assume we are given $ n$ nonlinear polynomial equations with $ l$ variables in the power basis, where $ n\geq l$ , and a box $ B = [a_{1}, b_{1}] \times [a_{2}, b_{2}] \times \ldots \times [a_{l}, b_{l}]$ , in which we need to determine the roots of the given system. In this case we first scale the box by performing an appropriate affine parameter transformation described above to the functions $ f_k$ , so that the box becomes $ [0,1]^l$ . Next we express the transformed nonlinear polynomial equations in the multivariate Bernstein basis using (4.18). Now we summarize the PP algorithm.

  1. Using the convex hull property, find a sub-box of $ [0,1]^l$ which contains all the roots. The essential idea behind the box generation scheme in this algorithm is to transform a complicated $ l+1$ -dimensional problem into a series of $ l$ two-dimensional problems. Suppose $ {\bf R}^{l+1}$ can be coordinatized with the $ u_1, u_2, \ldots, u_{l+1}$ axes; we can then employ these steps:

    1. Project the $ {\bf v}_I^{(k)}$ of all of the $ {\bf f}_k$ into $ l$ different coordinate planes; specifically, the $ (u_1, u_{l+1})$ -plane, the $ (u_2, u_{l+1})$ -plane, and so on, up to the $ (u_l, u_{l+1})$ plane.

    2. In each one of these planes,

      1. Construct $ n$ two-dimensional convex hulls. The first is the convex hull of the projected control points of $ {\bf f}_1$ , the second is from $ {\bf f}_2$ and so on.

      2. Intersect each convex hull with the horizontal axis (that is, $ u_{l+1} = 0$ ). Because the polygon is convex, the intersection may be either a closed interval (which may degenerate to a point) or empty. If it is empty, then no root of the system exists within the given search box.

      3. Intersect the intervals with one another. Again, if the result is empty, no root exists within the given search box.

    3. Construct an $ l$ -dimensional box by taking the Cartesian product of each one of these intervals in order. In other words, the $ u_1$ side of the box is the interval resulting from the intersection in the $ (u_1, u_{l+1})$ -plane, and so forth.

  2. Using the scaling relationship between our current box and the initial box of search, see if the new sub-box represents a sufficiently small box in $ {\bf R}^l$ . If it does not, then go to step 3. If it does, then check the convex hulls of the hypersurface in the new box. If the convex hulls cross each variable axis, conclude that there is a root or at least an approximate root in the new box, and put the new box into a root list. Otherwise the new box is discarded.

  3. If any dimension of this sub-box is not much smaller than 1 unit in length (i.e., the box has not decreased much in size along one or more sides), split the box evenly along each dimension which is causing trouble (not reducing in size). Continue on to the next iteration with several independent sub-problems.

  4. If none of the boxes is left, then the root-finding process is over. Otherwise, perform an appropriate affine parameter transformation to the functions $ f_k$ , so that the box becomes $ [0,1]^l$ , and go back to step 1 for each new box. This transformation can be performed with the multivariate de Casteljau subdivision algorithm which is an extension of similar algorithms for 1 and 2 dimensions given in [92]. However, keep track of the scaling relationship between this box and the initial box of search.

If we assume that each equation in (4.14) is of degree $ m$ in each variable and the system is $ n$ -dimensional, then the total asymptotic time per step is of $ O(nlm^{l+1})$ . The number of steps depends primarily on the accuracy required [392]. The Projected Polyhedron algorithm achieves quadratic convergence in one dimension, while for higher dimensions, it exhibits at best linear convergence [392]. Once roots have been isolated via the PP algorithm, local quadratically convergent Newton-type algorithms can be used to compute the roots to high precision more efficiently. An extension of the algorithm described above for a set of simultaneous piecewise polynomial nonlinear equations expressed in terms of tensor product B-splines can be found in [132]. A novel feature of this extension is the normalization of the equations in the range [-1,1] and normalization of the knot vector in each subdomain in range [0,1] at each iteration step of the process to capitalize on the higher density of floating point numbers in this range, thereby improving numerical robustness of the algorithm.

Because the PP algorithm depends only on the convex hull property and ability to perform subdivision and multiplication, in theory one could implement the algorithm for rational B-spline entities without subdividing them into their rational Bézier components. Subdivision algorithms for B-splines are well-known, and Mørken [270] has developed an algorithm for multiplying two piecewise polynomials expressed in the B-spline basis. However, Zhou et al. [461] indicate that this approach sometimes tends to be more time-consuming than subdividing into rational polynomials and applying direct algebraic operations of addition and multiplication of two Bernstein forms (see Sect. 1.3.2). Piegl and Tiller [315] provide detailed description of the procedures that can handle algebraic operators of NURBS curves and surfaces such as dot and cross products, sum/difference and derivative operators. They start with decomposing the B-splines into their Bézier components using knot insertion, and applying the algebraic operators to the Bézier functions and finally recomposing the resulting Bézier functions into B-spline form using knot removal.

next up previous contents index
Next: 4.5 Auxiliary variable method Up: 4. Nonlinear Polynomial Solvers Previous: 4.3.3 Subdivision Methods   Contents   Index
December 2009