Computer Aided Geometric Design 25 (2008) 529–541 www.elsevier.com/locate/cagd
A unified approach to resultant matrices for Bernstein basis polynomials Joab R. Winkler Department of Computer Science, The University of Sheffield, Regent Court, 211 Portobello Street, Sheffield S1 4DP, United Kingdom Received 22 November 2006; received in revised form 26 September 2007; accepted 27 September 2007 Available online 2 October 2007
Abstract Resultant matrices can be used to compute the intersection points of curves that are defined by polynomials. They were originally developed for polynomials expressed in the power (monomial) basis, but the recent development of resultant matrices for Bernstein basis polynomials has increased their use in computer aided geometric design, for which the Bernstein basis is the standard polynomial basis. In this paper, the equations that relate the Sylvester, Bézout and companion resultant matrices for Bernstein basis polynomials are derived, thereby establishing their equivalence. It is shown that these equations are more complicated than their power basis equivalents. © 2007 Elsevier B.V. All rights reserved. Keywords: Resultant matrices; Bernstein basis polynomials
1. Introduction The resultant of two polynomials is equal to the determinant of a resultant matrix, the entries of which are functions of their coefficients, and this can be used to compute the intersection points of curves (Bini et al., 2005; Goldman et al., 1984; De Montaudouin and Tiller, 1984). The recent development of resultant matrices for Bernstein basis polynomials has increased their use in computer aided geometric design because a potentially ill-conditioned transformation of polynomials from the Bernstein basis to the power basis is not required (Bini and Gemignani, 2004; Winkler, 2003; Winkler and Goldman, 2003). There exist several different resultant matrices, and they may be considered equivalent because their rank deficiency is equal to the degree of the greatest common divisor (GCD) of the polynomials, and the coefficients of the GCD can be obtained by performing an LU or QR decomposition on the resultant matrix. In this paper, the Bézout, companion and Sylvester resultant matrices for Bernstein basis polynomials are considered and the equations that relate them are derived, thereby establishing their equivalence. Section 2 contains a review of the Sylvester, Bézout and companion resultant matrices for power and Bernstein basis polynomials, and an example is included in order to show the form of each of the matrices. The equations that define the transformations of these matrices between their power and Bernstein forms are considered in Section 3, and these expressions are used in Section 4 to develop the formulae that relate these resultant matrices for Bernstein E-mail address:
[email protected]. 0167-8396/$ – see front matter © 2007 Elsevier B.V. All rights reserved. doi:10.1016/j.cagd.2007.09.004
530
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
basis polynomials. Section 5 considers some computational properties of resultant matrices, and Section 6 contains a summary of the paper. 2. Resultant matrices This section reviews the Sylvester, Bézout and companion resultant matrices for power and Bernstein basis polynomials. Let f (x) and g(x) be polynomials that are expressed as linear combinations of the Bernstein basis functions,1 and let p(x) and q(x) be their power basis forms respectively, m m m f (x) = ai ci x m−i , (1 − x)m−i x i and p(x) = i i=0
and g(x) =
i=0
n n (1 − x)n−i x i bi i
and q(x) =
i=0
n
di x n−i ,
i=0
where f (x) ≡ p(x) and g(x) ≡ q(x). 2.1. Bernstein basis resultant matrices The Sylvester resultant matrix, which is of order (m + n) × (m + n), of the polynomials f (x) and g(x) is (Winkler, 2005; Winkler and Goldman, 2003), ⎡ a m a m a m ⎤ · 0 0 1 1 2 2 ⎢ a0 m0 a1 m1 a2 m2 · ⎥ ⎢ ⎥ m m ⎢ ⎥ ⎢ a0 0 a1 1 · ·⎥ ⎢ ⎥ ⎢ · · ·⎥ ⎥ D, (1) S (b) (f, g) = ⎢ ⎢b n ⎥ b1 n1 b2 n2 · ⎢ 0 0 ⎥ ⎢ ⎥ ⎢ b0 n0 b1 n1 b2 n2 · ⎥ ⎢ ⎥ ⎣ ⎦ b0 n0 b1 n1 · · · · · where
1 1 1 1 (m+n)×(m+n) · · · D = diag (m+n−1) (m+n−1) . (m+n−1) (m+n−1) ∈ R 0
1
m+n−2
m+n−1
The Bézout resultant matrix Z (b) (f, g) of the polynomials f (x) and g(x) is of order r × r, where r = max(m, n). It is assumed that m n, and thus the polynomial g(x) is degree elevated (m − n) times (Winkler, 2006). The elements (b) zij , i, j = 1, . . . , m, of Z (b) (f, g) satisfy m f (x)g(y) − f (y)g(x) m−1 (m−1)−(i−1) i−1 m − 1 (1 − x) (1 − y)(m−1)−(j −1) y j −1 , zij x = x−y i−1 j −1 i,j =1
and expressions for them have been derived by Bini and Gemignani (2004), m (b) zi,1 = (ai b0 − a0 bi ), 1 i m, i m2 j (m − i) (b) zi,j +1 = (ai bj − aj bi ) + zi+1,j , 1 i, j m − 1, i(m − j ) i(m − j ) m (b) zm,j +1 = (am bj − aj bm ), 1 j m − 1. (m − j ) 1 For simplicity, f (x) and g(x) will be referred to in the sequel as Bernstein basis polynomials.
(2)
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
531
The eigenvalues of a square matrix A are the solutions of |A − λI | = 0, and the expansion of the determinant yields a polynomial equation in λ whose coefficients are non-trivial functions of the elements of A, such that their exact form is not obvious. By contrast, for a given polynomial whose coefficients and degree are specified, there exists an infinite number of matrices whose eigenvalues are equal to the roots of the polynomial. There is, however, a matrix whose eigenvalues are identically equal to the roots of the given polynomial, and a closed form expression for each of its elements, in terms of the coefficients of the polynomial, exists. This matrix C is called the companion matrix of a polynomial, and its form for the power basis polynomial p(λ) is specified in (5). It is readily verified that p(λ) = |C − λI |, and that the form of C is immediately obvious from the coefficients of p(λ). The polynomial matrix q(C) ∈ Rm×m is a companion resultant matrix because its determinant is identically zero if and only if the degree of the GCD of p(x) and q(x) is greater than zero. It is noted that C is the companion matrix of p(x), and that the coefficient of C n−i is di , which is the ith coefficient of q(x). More generally, the rank loss of q(C) is equal to the degree of the GCD of p(x) and q(x), and the coefficients of the GCD can be obtained by reducing q(C) to upper triangular form (Barnett, 1983). The discussion in the previous paragraphs considered the power basis polynomials p(x) and q(x), and a companion resultant matrix can also be developed for Bernstein basis polynomials (Winkler, 2003). Specifically, a companion resultant matrix for the polynomials f (x) and g(x) is g(M) ∈ Rm×m where M ∈ Rm×m is the companion matrix of f (x), and a closed form expression for its elements, in terms of the coefficients of g(x), exists (Winkler, 2003). The matrix g(M) possesses exactly the same properties as its power basis equivalent because its rank loss is equal to the degree of the GCD of f (x) and g(x), and the coefficients of their GCD can be obtained by reducing g(M) to upper triangular form. Example 2.1. Consider the Bernstein polynomials f (x) and g(x), 5 3 1 3 3 3 3 (1 − x)2 x − (1 − x)x 2 + x , f (x) = 3 (1 − x)3 − 6 1 2 2 3 0 and g(x) = 2
2 2 2 3 2 (1 − x)x + x , (1 − x)2 − 2 1 2 0
whose GCD is g(x) because 1 1 f (x) = − g(x) −3(1 − x) − 2x = − g(x)(x − 3). 2 2 It follows from (1) that the Sylvester resultant matrix S (b) (f, g) is ⎡ ⎤ ⎡3 −5 ⎡3 −5 −3 ⎤ 1 0 0 0 0 8 1 0 2 2 3 ⎢0 1 0 0 0⎥ ⎢0 ⎢0 3 −5 −3 1⎥⎢ ⎥ ⎢ 4 4 ⎢ ⎥ ⎢ ⎥⎢ 2 2 (b) 1 ⎢ ⎢ ⎥ = ⎢2 −3 ⎥ S (f, g) = ⎢ 2 −3 0 0 0 0 ⎥ ⎢ 6 4 1 0 0⎥⎢ ⎥ ⎢ ⎣ ⎦⎢ 1 0 2 −3 1 0 ⎣ 0 0 0 14 0 ⎦ ⎣ 0 2 0 0 2 −3 1 0 0 0 0 1 0 0 The calculation of the Bézout three, 3 g(x) = 2 (1 − x)3 − 0
3 4
1 4
5 − 12
− 38
1 6
0
− 12 1 3
1 4
− 34
0⎤ 1⎥ ⎥ ⎥ 0⎥ ⎥. ⎥ 0⎦ 1
resultant matrix (2) requires that g(x) be degree elevated to a polynomial of degree 1 3 2 3 3 3 2 2 (1 − x) x − (1 − x)x + x , 3 1 3 2 3
and thus the Bézout resultant matrix of f (x) and g(x) is ⎤ ⎡ 3 −2 −1 2 ⎢ 3 ⎥ Z (b) (f, g) = ⎣ 32 − 98 4 ⎦. −1
− 14
− 12
The companion matrix M of f (x) is (Winkler, 2003),
532
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
⎡
1
5 − 18
3
11 6 − 11 6
⎢ M = ⎣ −3
1 9
⎤
⎥ − 13 ⎦ , 4 3
and thus the companion resultant matrix of f (x) and g(x) is
⎡ 4 ⎤ 8 −4 3 2 2 3 2 g(M) = 2 (I − M)2 − (I − M)M + M 2 = ⎣ −36 18 −6 ⎦ . 0 2 1 2 54 −27 9
2.2. Power basis resultant matrices The Sylvester resultant matrix, which is of order (m + n) × (m + n), of the polynomials p(x) and q(x) is ⎡ ⎤ c0 c1 c2 · c0 c1 c2 · ⎥ ⎢ ⎢ ⎥ c0 c1 · · ⎥ ⎢ ⎢ ⎥ · · ·⎥ ⎢ S (p) (p, q) = ⎢ ⎥. ⎢ d0 d1 d2 · ⎥ ⎢ ⎥ d0 d1 d2 · ⎥ ⎢ ⎣ ⎦ d0 d1 · · · · ·
(3)
The Bézout resultant matrix Z (p) (p, q) of p(x) and q(x) is of order m × m since it is assumed that m n, and it (p) is therefore necessary to pad g(x) by (m − n) leading zeros. The elements zij , i, j = 1, . . . , m, of Z (p) (p, q) are (Barnett, 1983, page 45), (p)
zi,j =
s
|cm−(i+j −k−1) dm−k |,
s = min(i − 1, j − 1),
i, j = 1, . . . , m,
(4)
k=0
where |ci dj | = ci dj − cj di . The companion resultant matrix of the polynomials p(x) and q(x) is q(C), where C is the companion matrix of p(x), ⎡ 0 1 0 ··· 0 ⎤ ⎢ ⎢ C=⎢ ⎢ ⎣
0 .. .
0 .. .
1 .. .
0 − ccm0
0 − cm−1 c0
0 − cm−2 c0
··· .. .
0 ⎥ .. ⎥ . ⎥ ⎥. ··· 1 ⎦ · · · − cc10
(5)
The matrix C and polynomial matrix q(C) are each of order m × m. Example 2.2. Consider the power basis forms p(x) and q(x), respectively, of the Bernstein polynomials f (x) and g(x) in Example 2.1, p(x) = −3x 3 +
25 2 23 x − x + 3, 2 2
and q(x) = 6x 2 − 7x + 2,
whose GCD is q(x). It follows from (3) and (4) that the Sylvester and Bézout resultant matrices of p(x) and q(x) are ⎡ −3 25 − 23 3 0⎤ 2 2 ⎡ ⎤ 25 −2 7 −6 ⎢ 0 −3 − 23 3⎥ ⎢ ⎥ 2 2 (p) 49 ⎥ ⎣ S (p) (p, q) = ⎢ 21 ⎦ , ⎢ 6 −7 2 0 0 ⎥ and Z (p, q) = 7 − 2 ⎣ ⎦ −6 21 −18 0 6 −7 2 0 0 0 6 −7 2
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
533
respectively, and the power basis companion resultant matrix q(C) is 2 −7 6 q(C) = 6 −21 18 . 18 −63 54 It is readily checked that the rank of each of these resultant matrices is 2, and that reduction to row echelon form yields the coefficients of q(x). 3. Basis transformation equations of resultant matrices This section reviews the transformation equations between the power and Bernstein forms of the Sylvester, Bézout and companion resultant matrices. The Sylvester resultant matrix is considered initially, after which the Bézout and companion resultant matrices are considered. It is shown in Winkler and Goldman (2003) that the Sylvester resultant matrices S (b) (f, g) and S (p) (p, q) are related by2 −1 −1 (p) S (b) (f, g) = Em,n Nm,n S (p, q)Jm+n Tm+n−1 ,
where
Fn−1 0nm ∈ R(m+n)×(m+n) , 0mn Fm−1
d d d Fd = diag ··· ∈ R(d+1)×(d+1) , 0 1 d
(6)
Em,n =
(7) (8)
the subscripts on the zero matrices indicate their order, Jn ∈ Rn×n is the reverse unit matrix, Tm ∈ R(m+1)×(m+1) is the transformation matrix between the power and Bernstein basis functions for polynomials of degree m, ⎡ m ⎤ m ⎡ ⎤ 0 (1 − x) 1 ⎢ m (1 − x)m−1 x ⎥ ⎢ 1 ⎥ ⎢ x ⎥ ⎢ ⎥ ⎢ . ⎥ . ⎢ ⎥ = ⎢ .. ⎥ , .. (9) Tm ⎢ ⎥ ⎢ ⎥ ⎢ m ⎥ ⎣ x m−1 ⎦ m−1 ⎣ m−1 (1 − x)x ⎦ m m xm x m and
Nm,n =
Jn Tn−1 0mn
0nm ∈ R(m+n)×(m+n) . Jm Tm−1
(10)
The elements ti,j , i, j = 1, . . . , m + 1, of Tm are ⎧ if i > j, ⎨0 −1 ti,j = (ji−1 ) ⎩ m if i j, (i−1) −1 and the elements ti,j , i, j = 1, . . . , m + 1, of Tm−1 are 0 if i > j, −1 m j −1 = ti,j (j −i) if i j. (−1) j −1 i−1
The transformation between the Bézout resultant matrix for power and Bernstein basis polynomials, Z (p) (p, q) and Z (b) (f, g) respectively, is (Bini and Gemignani, 2004), 2 Eq. (6) contains the reverse unit matrix J m+n , which is not present in (11) in Winkler and Goldman (2003). This difference arises because the definitions of the power basis polynomials differ slightly between this work and Winkler and Goldman (2003).
534
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541 T Z (b) (f, g) = Tm−1 Z (p) (p, q)Tm−1 ,
(11)
where each matrix is of order m × m and Tm is defined in (9). The transformation of the companion resultant matrix between the power and Bernstein bases is (Winkler, 2004), q(C) = Hg(M)H −1 ,
(12)
where the elements hij , i, j = 1, . . . , m, of the upper triangular matrix H are ⎧ if i > j, ⎨0 j −1 hij = m+1−j ( i−1 ) ⎩ ( m ) m−1 if i j, ( i−1 ) −1 are and the elements h−1 ij , i, j = 1, . . . , m, of H
h−1 ij =
0 m m−i (−1)(j −i) i−1 j −i
if i > j, if i j.
Example 3.1. Consider the polynomials in Examples 2.1 and 2.2, for which m = 3 and n = 2. Consider initially the transformation (6) between the power and Bernstein forms of the Sylvester resultant matrix, ⎡ ⎤ 0 0 0 0 1 ⎢0 0 0 1 1⎥
4 ⎢ ⎥ 1 ⎢ ⎥ −1 E3,2 J5 T4 = ⎢ 0 0 16 21 1 ⎥ , = diag 1 1 1 1 , ⎢ ⎥ 2 ⎣0 1 1 3 1⎦ 4 2 4 1 1 1 1 1 and
⎡
−1 N3,2 =
T1−1 J2 0
−1 ⎢ 1 0 ⎢ =⎢ 0 −1 ⎣ T2 J 3 0 0
⎤ 1 0 0 0 0 0 0 0⎥ ⎥ 0 1 −2 1 ⎥ . ⎦ 0 −2 2 0 0 1 0 0
It is readily checked that the right-hand side of (6) is equal to S (b) (f, g), as required. Consider now the transformation of the Bézout resultant matrix between the power and Bernstein bases. Since ⎡ ⎤ 1 1 1 1 −2 1 T2 = ⎣ 0 12 1 ⎦ and T2−1 = 0 2 −2 , 0 0 1 0 0 1 and, from Example 2.2, ⎡ −2 7 Z (p) (p, q) = ⎣ 7 − 49 2 −6 21 it follows that (11) yields ⎡ 1 0 Z (b) (f, g) = ⎣ 1 12 1 1
⎤ −6 21 ⎦ , −18
⎤⎡ 0 −2 7 0 ⎦ ⎣ 7 − 49 2 1 −6 21
⎤⎡ −6 1 21 ⎦ ⎣ 0 −18 0
1 1 2
0
⎤ ⎡ −2 3 1 2 ⎢ 1 ⎦ = ⎣ 32 − 98 3 1 −1 4
−1 3 4
⎤ ⎥ ⎦,
− 12
which agrees with the Bézout resultant matrix Z (b) (f, g) in Example 2.1. Finally, consider the transformation of the companion resultant matrix between the power and Bernstein bases (12). Since
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
⎡
1
⎢ H = ⎣0 0
2 3 1 3
0
1 3 1 3 1 3
⎤ ⎥ ⎦
535
and H
−1
1 −2 1 = 0 3 −3 , 0 0 3
it follows from (12) that
1 −2 1 H −1 q(C)H = 0 3 −3 0 0 3
⎡ 1 2 −7 6 ⎢ 6 −21 18 ⎣ 0 18 −63 54 0
2 3 1 3
0
1 3 1 3 1 3
⎤
⎡
4 ⎤ 8 −4 3 ⎥ ⎣ ⎦ = −36 18 −6 ⎦ , 54 −27 9
which agrees with the Bernstein basis companion resultant matrix in Example 2.1. 4. The equations that relate the resultant matrices This section considers the equations that relate the Sylvester, Bézout and companion resultant matrices in the Bernstein basis. The equivalent formulae for the power basis are developed in Barnett (1983), and it is shown that the formulae for the Bernstein basis are more complex. 4.1. The Bézout and companion resultant matrices It is shown in Barnett (1983), page 46, that the Bézout and companion resultant matrices in the power basis are related, Z (p) (p, q) = Rq(C),
(13)
∈ Rm×m , the companion matrix of
where C contain the coefficients of p(x), ⎡c c c m−1
m−2
m−3
⎢ cm−2 ⎢ ⎢ cm−3 ⎢ R=⎢ . ⎢ .. ⎢ ⎣ c 1 c0
cm−3 cm−4 .. . c0 0
cm−4 cm−5 .. . 0 0
... ... ... .. . ... ...
c1 c0 0 .. . 0 0
p(x), is defined in (5), and R is a Hankel matrix whose columns (or rows) c0 ⎤ 0⎥ ⎥ 0⎥ ⎥ . .. ⎥ . ⎥ ⎥ 0⎦ 0
Example 4.1. The Hankel matrix R for the polynomial p(x), which is defined in Example 2.2, is ⎡ 23 25 ⎤ −2 −3 2 R = ⎣ 25 −3 0 ⎦ , 2 −3 0 0 and thus ⎡ 23 25 ⎤ ⎤ ⎡ −2 −2 −3 2 7 −6 −7 6 2 Rq(C) = ⎣ 25 21 ⎦ , −3 0 ⎦ 6 −21 18 = ⎣ 7 − 49 2 2 18 −63 54 −6 21 −18 −3 0 0 which is equal to Z (p) (p, q), as required. Eq. (13) shows that the Bézout and companion resultant matrices in the power basis are related by the Hankel matrix R, and since R contains the coefficients of p(x), it is necessary to determine its form for Bernstein basis (p) polynomials. In particular, if r(x) is the constant polynomial r(x) = 1, then the elements zi,j of Z (p) (p, r) are calculated from (4), (p)
zi,j =
s k=0
|cm−(i+j −k−1) rm−k | = cm−(i+j −1) ,
s = min(i − 1, j − 1),
536
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
because the only non-zero coefficient of r(x) is rm = 1. Since this is the expression for element (i, j ) of R, it follows that Z (p) (p, r) = R. The interpretation of R as a Bézout resultant matrix in the power basis allows (11) to be used to determine its Bernstein basis form. In particular, it follows from this equation that, with g(x) = q(x) = r(x) = 1, T T Z (b) (f, r) = Tm−1 Z (p) (p, r)Tm−1 = Tm−1 RTm−1 ,
and thus −T −1 R = Tm−1 Z (b) (f, r)Tm−1 .
(14)
It follows from (12), (13) and (14) that −T −1 Z (p) (p, q) = Rq(C) = Tm−1 Z (b) (f, r)Tm−1 Hg(M)H −1 ,
and thus −1 T Tm−1 Z (p) (p, q)Tm−1 = Z (b) (f, r) Tm−1 H g(M) H −1 Tm−1 , and hence from (11), −1 Z (b) (f, g) = Z (b) (f, r) Tm−1 H g(M) H −1 Tm−1 .
(15)
This equation expresses the Bézout resultant matrix Z (b) (f, g) as a function of the Bézout resultant matrix of f (x) and the constant polynomial r(x) = 1, and the companion resultant matrix g(M), and it therefore follows that this equation has a more complex form than its power basis equivalent (13). Eq. (15) can be simplified because T , H, T −1 and H −1 are upper triangular matrices for which closed form expressions for their elements are known. Eq. (9) shows that Tm−1 is a basis transformation matrix, and it is shown in Winkler (2004) that H is also a basis transformation matrix, ⎡ m (1 − x)m−1 ⎤ ⎡ ⎤ 1 0 ⎢ m (1 − x)m−2 x ⎥ ⎥ ⎢ x ⎥ ⎢ 1 ⎥ ⎢ . ⎥ ⎢ . ⎥ = ⎢ .. ⎥ , ⎢ H⎢ .. ⎥ ⎢ ⎥ ⎥ ⎣ x m−2 ⎦ ⎢ m ⎣ m−2 (1 − x)x m−2 ⎦ m m−1 x m−1 m−1 x and it therefore follows that on changing m into (m − 1) in (9), ⎡ m−1 ⎡ m m−1 ⎤ m−1 ⎤ 0 (1 − x) 0 (1 − x) ⎢ m−1 ⎥ ⎢ m (1 − x)m−2 x ⎥ ⎢ 1 (1 − x)m−2 x ⎥ ⎥ ⎢ 1 ⎢ ⎥ ⎥ ⎢ ⎢ ⎥ −1 .. .. ⎥. ⎢ = T H ⎢ ⎥ m−1 ⎥ ⎢ . . ⎢ ⎥ ⎥ ⎢ m ⎢ m−1 (1 − x)x m−2 ⎥ m−2 ⎦ ⎣ (1 − x)x ⎣ m−2 ⎦ m−2 m−1 m−1 m m−1 m−1 x m−1 x −1 This equation shows that Tm−1 H is a diagonal matrix G,
m+1−i m m−1 −1 H= = diag 1 G = Tm−1 m m i=1
...
2 m
1 , m
(16)
and thus (15) can be written as Z (b) (f, g) = Z (b) (f, r)Gg(M)G−1 .
(17)
This equation relates the Bézout resultant matrix of the Bernstein basis polynomials f (x) and g(x), and the companion resultant matrix g(M), where M is the companion matrix of f (x). It involves the Bézout resultant matrix of f (x) and the constant polynomial r(x) = 1.
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
537
Example 4.2. The determination of the Bézout resultant matrix Z (b) (f, r) requires that r(x) be degree elevated to a polynomial of degree three, and it is clear that the coefficients of the degree elevated form of r(x) are ri = 1, i = 0, . . . , 3. It therefore follows from (2) and (16) that ⎤ ⎡ 23 ⎡ ⎤ − 2 − 21 −2 1 0 0 4 ⎥ ⎢ 2 1 11 ⎣0 3 0 ⎦, Z (b) (f, r) = ⎣ − 21 4 4 4 ⎦ and G = 9 11 0 0 13 −2 4
2
respectively, and thus (17) yields ⎤ ⎡ 23 −2 ⎡ 1 − 2 − 21 4 ⎢ 1 11 ⎥ ⎣ 0 Z (b) (f, g) = ⎣ − 21 4 4 4 ⎦ 9 11 0 −2 4 2 ⎡ ⎤ 3 −1 −2 2 ⎢ 3 3 ⎥ 9 = ⎣ 2 −8 4 ⎦, −1
3 4
0 2 3
0
⎤⎡ 4 ⎤⎡1 0 8 −4 3 0 ⎦ ⎣ −36 18 −6 ⎦ ⎣ 0 1 54 −27 9 0 3
0 3 2
0
⎤ 0 0⎦ 3
− 12
which agrees with the result in Example 2.1. 4.2. The Sylvester and companion resultant matrices It is shown in Barnett (1983), page 32, that the Sylvester and companion resultant matrices in the power basis are related by
In 0nm (p) Pnm Pnn S (p, q) = , (18) K Im 0mn Jm q(C)Jm where K ∈ Rm×n is an upper triangular Toeplitz matrix, and the Sylvester resultant matrix is partitioned as
Pnn Pnm , S (p) (p, q) = Pmn Pmm
(19)
where Pab ∈ Ra×b . It follows from (18) that K satisfies KPnn + Pmn = 0 and KPnm + Pmm = Jm q(C)Jm .
(20)
Eq. (18) relates the power basis forms of the Sylvester and companion resultant matrices, and it is desired to determine the equivalent formula for the Bernstein basis. It follows from (6) that −1 S (p) (p, q) = Nm,n Em,n S (b) (f, g)Tm+n−1 Jm+n ,
where, from (7) and (10),
Jn Tn−1 Fn−1 Nm,n Em,n = 0mn
(21)
0nm . Jm Tm−1 Fm−1
−1 If the product Tm+n−1 Jm+n is partitioned as −1
−1 −1 0 T Tnm Tnn Im mn −1 = nn Tm+n−1 Jm+n = −1 −1 −1 Tmm Tmn In 0nm Tmn −1 where Tab ∈ Ra×b , and S (b) (f, g) is partitioned as
Bnn Bnm (b) , S (f, g) = Bmn Bmm
where Bab ∈ Ra×b , it follows from (19) and (21) that
−1 Tnm −1 , Tmm
538
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
−1 −1 , Pnn = Jn Tn−1 Fn−1 Bnn Tnn + Bnm Tmn −1 −1 , + Bnm Tmm Pnm = Jn Tn−1 Fn−1 Bnn Tnm −1 −1 Pmn = Jm Tm−1 Fm−1 Bmn Tnn + Bmm Tmn , −1 −1 . + Bmm Tmm Pmm = Jm Tm−1 Fm−1 Bmn Tnm The substitution of these expressions into (20) yields −1 K = −Pmn Pnn
−1 −1 −1 −1 −1 −1 −1 Bnn Tnn = −Jm Tm−1 Fm−1 Bmn Tnn + Bmm Tmn + Bnm Tmn Fn−1 Tn−1 Jn ,
and it follows from (12) and (20) that KPnm + Pmm = Jm Tm−1 Fm−1 −1 −1 −1 −1 − Bmn Tnn × Bmn Tnm + Bmm Tmm + Bmm Tmn −1 −1 −1 −1 −1 × Bnn Tnn Bnn Tnm + Bnm Tmn + Bnm Tmm = Jm Hg(M)H −1 Jm . Hence
−1 −1 −1 −1 g(M) = H −1 Tm−1 Fm−1 Bmn Tnm − Bmn Tnn + Bmm Tmm + Bmm Tmn −1 −1 −1 −1 −1 × Bnn Tnn Bnn Tnm Jm H, + Bnm Tmn + Bnm Tmm
and this expression can be simplified by defining the matrix W
Wnn Wnm W := Wmn Wmm −1
−1 Tnn Tnm Bnn Bnm = −1 T −1 Bmn Bmm Tmn mm
−1 −1 −1 + B T −1 Bnn Tnm Bnn Tnn + Bnm Tmn nm mm = −1 −1 −1 −1 , Bmn Tnn + Bmm Tmn Bmn Tnm + Bmm Tmm
∈ R(m+n)×(m+n)
(22) as
(23)
where Wab ∈ Ra×b . Eq. (22) therefore yields, using (16), −1 g(M) = H −1 Tm−1 Fm−1 Wmm − Wmn Wnn Wnm Jm H = H −1 Tm−1 Fm−1 (W/Wnn )Jm H = G−1 Fm−1 (W/Wnn )Jm H = Lm−1 (W/Wnn )Jm H,
(24)
where (W/Wnn ) is the Schur complement of Wnn in W , and Lm−1 ∈ Rm×m is the diagonal matrix whose elements are calculated from (8) and (16), m
m m m m = diag ··· . Lm−1 = G−1 Fm−1 = i − 1 i=1 0 1 m−1 Eq. (24) relates the Sylvester and companion resultant matrices for Bernstein basis polynomials, where the Sylvester resultant matrix is defined by its blocks Bnn , Bnm , Bmn and Bmm . The equation can be written in a slightly different form,
Wnn Wnm Wnn 0nm Wnm In = , (25) −1 −1 J −Wmn Wnn Im Wmn Wmm 0mn L−1 m m−1 g(M)H where, from (23),
−1 −1 Wnn Wnm Tnn Tnm (b) = S (f, g) −1 T −1 . Wmn Wmm Tmn mm
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
539
Example 4.3. The Schur complement (W/W22 ) for the polynomials f (x) and g(x) is 24 −28 8 (W/W22 ) = −36 42 −12 , 54 −63 18 and the matrices L2 and J3 H are, respectively, ⎡ 0 0 ⎢ L2 = diag[1 3 3] and J3 H = ⎣ 0 13 1
2 3
1 3 1 3 1 3
⎤ ⎥ ⎦.
It is readily checked that the right-hand side of (24) yields g(M), which is specified in Example 2.1. 4.3. The Bézout and Sylvester resultant matrices This section establishes the equation that relates the Bézout and Sylvester resultant matrices for Bernstein basis polynomials. Its derivation follows easily from (24) because Gg(M)G−1 = Fm−1 (W/Wnn )Jm H G−1 , and thus from (17), (b) −1 Z (f, r) Z (b) (f, g) = Fm−1 (W/Wnn )Jm H G−1 . It therefore follows from (16) that Z (b) (f, g) = Z (b) (f, r)Fm−1 (W/Wnn )Jm Tm−1 ,
(26)
which establishes the relation between the Bézout and Sylvester resultant matrices for Bernstein basis polynomials. This equation can be written as
−1
−1 0nm (b) Tnn Tnm In S (f, g) −1 −1 T −1 −Wmn Wnn Im Tmn mm
Wnm Wnn . (27) = −1 −1 0mn Fm−1 (Z (b) (f, r))−1 Z (b) (f, g)Tm−1 Jm Example 4.4. The Bézout resultant matrix Z (b) (f, r) and the Schur complement (W/W22 ) for the polynomials f (x) and g(x) are specified in Examples 4.2 and 4.3, respectively. The matrices F2 and J3 T2 are respectively, ⎡ ⎤ 0 0 1 F2 = diag[1 2 1] and J3 T2 = ⎣ 0 12 1 ⎦ , 1 1 1 and it is readily checked that the right-hand side of (26) yields Z (b) (f, g), which is specified in Example 2.1. 5. Computational issues Section 4 considered the equations that relate the Sylvester, Bézout and companion resultant matrices for Bernstein basis polynomials, and this section addresses some computational issues of resultant matrices. It is well known that the Bernstein basis is numerically superior to the power basis in the interval [0, . . . , 1], and that a basis transformation is not recommended because it is ill-conditioned. This ill-conditioned property is more severe for the transformation of a resultant matrix between these bases. In particular, it is shown in Winkler (2005) that the condition number of the transformation between the power and Bernstein forms of the Bézout resultant matrix, which is stated in (11), is equal to the square of the condition number of Tm−1 , and the equivalent condition number for the transformation (12) for the companion resultant matrix is equal to the square of the condition number of H .
540
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
Comparisons of the numerical properties of the power and Bernstein forms of resultant matrices have been performed and results that are consistent with the superior numerical properties of the Bernstein basis have been obtained. For example, Bini and Gemignani (2004), and Winkler (2003), compare the power and Bernstein forms of the Bézout and companion resultant matrices, respectively, for the computation of the GCD of two polynomials, and it is shown that the Bernstein form of the resultant matrix yields consistently better numerical answers than does its power basis equivalent. This comparison cannot be extended to the Sylvester resultant matrix because the coefficients of the polynomials are decoupled in this matrix. In particular, it is clear from (1) and (3) that S (b) (f, αg) = αS (b) (f, g)
and S (p) (p, αq) = αS (p) (p, q),
(28)
where α is an arbitrary non-zero scalar, and thus, for example, the condition number and distance to singularity, and the QR and LU decompositions, of the Sylvester resultant matrix for both bases can be changed by merely scaling one of the polynomials. This numerical property must be considered when computations with the Sylvester resultant matrix are performed. The property (28) of the Sylvester resultant matrix is not shared by the companion and Bézout resultant matrices. Specifically, it is readily verified that Z (b) (f, αg) = αZ (b) (f, g)
and Z (p) (p, αq) = αZ (p) (p, q),
and thus the Bézout resultant matrix is scale invariant because scaling one or both polynomials scales, for example, the singular values of Z (b) (f, g) and Z (p) (p, q), such that the ratio of any two singular values is constant. Since the companion matrices C and M require that the coefficients of x m in the power and Bernstein forms, respectively, of the polynomials be one, it necessarily follows that the companion resultant matrix is also scale invariant. The scale invariance of the companion and Bézout resultant matrices must be compared with the lack of scale invariance of the Sylvester resultant matrix (28). This difference arises because every entry in the Sylvester resultant matrix is linear in the coefficients of exactly one, and not both, polynomials, that is, the coefficients of the polynomials are decoupled in this matrix. Every entry in the Bézout resultant matrix is, however, exactly bilinear in the coefficients of the polynomials, with the consequence that scaling the coefficients of one polynomial scales the entire matrix. It therefore follows that the effect of α must be included in a study of the numerical properties of the Sylvester resultant matrix, but not in a study of the Bézout and companion resultant matrices. Resultant matrices are used in several applications, including fast polynomial arithmetic (Mourrain and Pan, 2000), polynomial root finding (Cardinal, 1996; Zeng, 2005), and the computation of the GCD of two polynomials (Barnett, 1983, 1988; Winkler and Allan, 2007a, 2007b), which arises in several applications, including computer aided geometric design, signal processing and control. Although there exists an extensive literature on the computation of the GCD of two polynomials, it is a non-trivial problem because a continuous change in the coefficients of one or both polynomials may cause a discontinuous change in the GCD, which implies that it is an ill-posed problem. It is well-known that even if the coefficients of the polynomials are specified exactly, roundoff errors are sufficient to cause totally incorrect answers. Furthermore, if the coefficients of the polynomials are inexact, which occurs frequently in practical applications, then even if the theoretically exact forms of the polynomials have a non-constant GCD, the given inexact polynomials are, with probability of almost 1, coprime. In this circumstance, an approximate GCD can be defined, that is, a non-constant GCD that is obtained by perturbing the given inexact polynomials by a small amount, such that the perturbed polynomials are not coprime. It is shown in Winkler and Allan (2007a, 2007b) that an approximate GCD can be calculated by computing perturbations of the Sylvester resultant matrix, such that its structure is preserved and a rank deficiency is introduced. Moreover, it is shown in these references that it is necessary to include the parameter α, whose role is defined in (28), in the computations in order to compute a valid approximate GCD. 6. Summary The equations that relate the Sylvester, Bézout and companion resultant matrices for Bernstein polynomials have been derived ((17), (25) and (27)), and examples have been given. These equations assume a more complicated form than their power basis equivalents because resultant matrices for Bernstein polynomials have a more complex structure than do resultant matrices for power basis polynomials. The equations between the Sylvester and companion resultant matrices, and the Sylvester and Bézout resultant matrices, (25) and (27) respectively, involve the Schur complement,
J.R. Winkler / Computer Aided Geometric Design 25 (2008) 529–541
541
and the equation between the Bézout and companion resultant matrices requires that a Hankel matrix be interpreted as a Bézout resultant matrix (17). A discussion of the applications of resultant matrices was included, and some computational issues were discussed. It was shown that it is necessary to introduce a scale parameter in the Sylvester resultant matrix in order to obtain good computational answers because the coefficients of the polynomials are decoupled in this resultant matrix, and thus one polynomial can be scaled relative to the other polynomial, with a consequent change in the computed solution. This feature of the Sylvester resultant matrix must be compared with the Bézout and companion resultant matrices, for which scaling is not required. References Barnett, S., 1983. Polynomials and Linear Control Systems. Marcel Dekker, New York. Barnett, S., 1988. Euclidean remainders for generalized polynomials. Linear Algebra Appl. 99, 111–122. Bini, D.A., Gemignani, L., 2004. Bernstein–Bezoutian matrices. Theoret. Comput. Sci. 315, 319–333. Bini, D.A., Gemignani, L., Winkler, J.R., 2005. Structured matrix methods for CAGD: An application to computing the resultant of polynomials in the Bernstein basis. Numer. Linear Algebra Appl. 12, 685–698. Cardinal, J.P., 1996. On two iterative methods for approximating the roots of a polynomial. In: Renegar, J., Shub, M., Smale, S. (Eds.), The Mathematics of Numerical Analysis, Lecture Notes in Applied Mathematics, vol. 32, pp. 165–188. Goldman, R.N., Sederberg, T.W., Anderson, D.C., 1984. Vector elimination: A technique for the implicitization, inversion and intersection of planar parametric rational polynomial curves. Comput. Aided Geom. Design 1, 327–356. De Montaudouin, Y., Tiller, W., 1984. The Cayley method in computer aided geometric design. Comput. Aided Geom. Design 1, 309–326. Mourrain, B., Pan, V.Y., 2000. Multivariate polynomials, duality and structured matrices. J. Complexity 16, 110–180. Winkler, J.R., 2003. A companion matrix resultant for Bernstein polynomials. Linear Algebra Appl. 362, 153–175. Winkler, J.R., 2004. The transformation of the companion matrix resultant between the power and Bernstein polynomial bases. Appl. Numer. Math. 48 (1), 113–126. Winkler, J.R., 2005. Numerical and algebraic properties of Bernstein basis resultant matrices. In: Dokken, T., Jüttler, B. (Eds.), Computational Methods for Algebraic Spline Surfaces. Springer-Verlag, Berlin, pp. 107–118. Winkler, J.R., 2006. The numerical condition of univariate and bivariate degree elevated Bernstein polynomials. J. Comp. Appl. Math. 191 (1), 32–49. Winkler, J.R., Allan, J.D., 2007a. Structured low rank approximations of the Sylvester resultant matrix for approximate GCDs of Bernstein polynomials. Submitted for publication. Winkler, J.R., Allan, J.D., 2007b. Structured total least norm and approximate GCDs of inexact polynomials. J. Comput. Appl. Math., in press. Winkler, J.R., Goldman, R.N., 2003. The Sylvester resultant matrix for Bernstein polynomials. In: Lyche, T., Mazure, M., Schumaker, L.L. (Eds.), Curve and Surface Design. Saint-Malo, 2002. Nashboro Press, Tennessee, pp. 407–416. Zeng, Z., 2005. Computing multiple roots of inexact polynomials. Math. Comput. 74, 869–903.