Structured matrix methods for polynomial computations Joab R. - - PowerPoint PPT Presentation

✬ ✫ ✩ ✪

Structured matrix methods for polynomial computations

Joab R. Winkler Department of Computer Science The University of Sheffield United Kingdom Structured Matrix Days, Limoges, France. May 2012

✬ ✫ ✩ ✪ Contents

• Introduction

– Difficulties of performing reliable computations on polynomials – The condition number of a root of a polynomial

• Multiple roots
• The geometry of multiple roots
• A simple polynomial root finder

– GCD computations – Polynomial division

• Examples
• Summary

1

✬ ✫ ✩ ✪

• 1. INTRODUCTION

Problems in which it is required to compute the greatest common divisor (GCD) of two polynomials, and the roots of polynomials:

• Computer-aided geometric design:

The calculation of the points of intersection of curves and surfaces

• Control theory
• Cryptography

2

✬ ✫ ✩ ✪

1.1 Difficulties of performing reliable computations on polynomials Commonly used methods and algorithms for computing the roots of a polynomial include:

• Bairstow, Graeffe, Jenkins-Traub, Laguerre, M¨

uller, Newton, QR decomposition Commonly used software for computing the roots of a polynomial include:

• MATHEMATICA, MATLAB, NAG, MPSOLVE

Satisfactory results are obtained if:

• The polynomial is of moderate degree
• The roots are simple and well-separated
• A good starting point in the iterative scheme is used

3

✬ ✫ ✩ ✪

Example 1.1 Consider the polynomial

y4 − 4y3 + 6y2 − 4y + 1 = (y − 1)4

whose root is y = 1 with multiplicity 4. MATLAB returns the roots 1.0002, 1.0000 ± 0.0002i, 0.9998

• Example 1.2 The roots of the polynomial (y − 1)100 were computed by MATLAB.

1 2 3 4 5 6 −4 −3 −2 −1 1 2 3 4 Real Imag

Figure 1.1: The computed roots of

(y − 1)100.

• 4

✬ ✫ ✩ ✪

Example 1.3 The computation of multiple roots in the presence of noise.

0.5 0.6 0.7 0.8 0.9 1 1.1 1.2 1.3 −0.08 −0.06 −0.04 −0.02 0.02 0.04 0.06 0.08

f(x) = (x−0.6)3(x−1)5 εc = 1e−007

Real Imag 0.7 0.8 0.9 1 1.1 −0.1 −0.05 0.05 0.1

f(x) = (x−0.7)3(x−1)5 εc = 1e−007

Real Imag 0.75 0.8 0.85 0.9 0.95 1 1.05 1.1 1.15 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15

f(x) = (x−0.8)3(x−1)5 εc = 1e−007

Real Imag 0.8 0.9 1 1.1 1.2 1.3 −0.2 −0.15 −0.1 −0.05 0.05 0.1 0.15

f(x) = (x−0.9)3(x−1)5 εc = 1e−007

Real Imag

Figure 1.2: The distribution of the roots of four perturbed polynomials.

• 5

✬ ✫ ✩ ✪

Example 1.4 Consider the polynomial (y − 1)100 and perturb the constant term by

2−100 g(y) = (y − 1)100 − 2−100

The roots of g(y) are

yk = 1 + 1 2

• cos

2kπ 100

• + i sin

2kπ 100

• ,

k = 0, . . . , 99

A change in the constant coefficient of 2−100 ≈ 10−32 causes a change in the magnitude of the roots from 1 to 1.5.

• 6

✬ ✫ ✩ ✪

• 2. MULTIPLE ROOTS

Consider the polynomial

f(y) =

m

• i=0

aiym−i

that has a root of multiplicity r

f (k)(y) = 0, f (0)(y) = f(y), k = 0, . . . , r − 1

• Why are multiple roots ill-conditioned?
• Are there conditions for which they are well-conditioned, and can they be utilised

in a stable algorithm for their computation?

7

✬ ✫ ✩ ✪

Example 2.1

εc =

componentwise noise amplitude componentwise signal amplitude = 10−8 exact exact computed computed relative root mult. root mult. error

• 6.7547e-01

4

• 6.7547000082e-01

4 1.2139725913e-09 5.7335e+00 6 5.7335000822e+00 6 1.4344694923e-08 2.1747e+00 7 2.1746999924e+00 7 3.5069931355e-09

• 9.5568e+00

10

• 9.5567996740e+00

10 3.4111255034e-08

• 6.5553e+00

11

• 6.5553001701e+00

11 2.5954947075e-08

• 8

✬ ✫ ✩ ✪

Example 2.2

εc =

componentwise noise amplitude componentwise signal amplitude = 10−8 exact exact computed computed relative root mult. root mult. error

• 1.1539e+00

4

• 1.1539000316e+00

4 2.7389266674e-08 4.0809e+00 5 4.0808998890e+00 5 2.7198581384e-08

• 2.1059e+00

6

• 2.1058999294e+00

6 3.2521823798e-08 3.6683e+00 7 3.6683000481e+00 7 1.3110114060e-08

• 9.6084e+00

13

• 9.6084001121e+00

13 1.1664214687e-08

• 9

✬ ✫ ✩ ✪

• 3. THE GEOMETRY OF MULTIPLE ROOTS

William Kahan (UC Berkeley) showed it is necessary to distinguish between a random perturbation, and a structured perturbation, applied to the coefficients of a polynomial

f(y) when considering the numerical condition of its roots.

Let the root yi of f(y) have multiplicity mi:

f(y) =

m

• i=0

aiym−i =

r

• i=1

(y − yi)mi,

r

• i=1

mi = m

Consider the perturbed polynomials g(y) and h(y):

g(y) =

m

• i=0

(ai + δai)ym−i, |δai| ≤ ε |ai| h(y) =

r

• i=1

(y − (yi + δyi))mi , |δyi| ≪ |yi|

10

✬ ✫ ✩ ✪

• Consider the polynomial g(y)

g(y) =

m

• i=0

(ai + δai)ym−i, |δai| ≤ ε |ai|

If the roots of g(y) are yi + δyi, then the roots yi are ill-conditioned.

∆yi ε ≫ 1, ∆yi = |δyi| |yi| , i = 1, . . . , m

• Consider the polynomial h(y)

h(y) =

r

• i=1

(y − (yi + δyi))mi , |δyi| ≪ |yi|

Then

∆yi ∆f = 0(1), ∆yi = |δyi| |yi| , i = 1, . . . , r, ∆f = ||f − h f

11

✬ ✫ ✩ ✪

Summary:

• If random perturbations are applied to the coefficients of a polynomial f(y), then

its multiple roots are ill-conditioned and break up into complex conjugate simple roots.

• If structured perturbations are applied to the coefficients of f(y), such that the

multiplicities of its roots are preserved, then its roots are stable. Kahan introduced the pejorative manifold of f(y):

f(y) =

r

• i=1

(y − yi)mi =

m

• i=0

aiym−i,

r

• i=1

mi = m

• The pejorative manifold of f(y) defines the constraints satisfied by the coeffi-

cients ai such that the multiplicities of its roots are mi.

12

✬ ✫ ✩ ✪

Example 3.1 Consider a quartic polynomial f(y) with real roots y1, y2, y3 and y4:

f(y) = y4 − (y1 + y2 + y3 + y4)y3 + (y1y2 + y1y3 + y1y4 + y2y3 + y2y4 + y3y4)y2 − (y1y2y3 + y1y2y4 + y1y3y4 + y2y3y4)y + y1y2y3y4

• If f(y) has a quartic root, then y1 = y2 = y3 = y4, then

f(y) = y4 − 4y1y3 + 6y2

1y2 − 4y3 1y + y4 1

The pejorative manifold of this quartic polynomial is

• −4y1

6y2

1

−4y3

1

y4

1

• 13

✬ ✫ ✩ ✪

• If f(y) has one treble root and one simple root, then y1 = y2 = y3 = y4, then

f(y) = y4 − (3y1 + y4)y3 + 3(y2

1 + y1y4)y2 − (y3 1 + 3y2 1y4)y + y3 1y4

The pejorative manifold of this quartic polynomial is

• −(3y1 + y4)

3(y2

1 + y1y4)

−(y3

1 + 3y2 1y4)

y3

1y4

• If f(y) has two double roots, then y1 = y2 and y3 = y4 = y1, then

f(y) = y4 − 2(y1 + y4)y3 + (y2

1 + 4y1y4 + y2 4)y2

−2(y2

1y4 + y1y2 4)y + y2 1y2 4

The pejorative manifold of this quartic polynomial is

• −2(y1 + y4)

(y2

1 + 4y1y4 + y2 4)

−2(y2

1y4 + y1y2 4)

y2

1y2 4

• 14

✬ ✫ ✩ ✪

• 4. A SIMPLE POLYNOMIAL ROOT FINDER

Gauss developed a polynomial root finder that is fundamentally different from the methods mentioned earlier. Let w1(y) be the product of all linear factors of f(y) Let w2(y) be the product of all quadratic factors of f(y)

· · ·

Let wi(y) be the product of all factors of degree i of f(y)

f(y) = w1(y)w2

2(y)w3 3(y) · · · wr r(y)

15

✬ ✫ ✩ ✪

Perform a sequence of GCD computations:

q1(y) =

GCD

• f(y), f (1)(y)
• =

w2(y)w2

3(y)w3 4(y) · · · wr−1 r

(y) q2(y) =

GCD

• q1(y), q(1)

1 (y)

• =

w3(y)w2

4(y)w3 5(y) · · · wr−2 r

(y) q3(y) =

GCD

• q2(y), q(1)

2 (y)

• =

w4(y)w2

5(y)w3 6(y) · · · wr−3 r

(y) q4(y) =

GCD

• q3(y), q(1)

3 (y)

• =

w5(y)w2

6(y)w3 7(y) · · · wr−4 r

(y)

. . . The sequence terminates at qr(y), which is a constant

16

✬ ✫ ✩ ✪

A set of polynomials hi(y), i = 1, . . . , r, is defined such that

h1(y) =

f(y) q1(y)

= w1(y)w2(y)w3(y) · · · h2(y) =

q1(y) q2(y)

= w2(y)w3(y) · · · h3(y) =

q2(y) q3(y)

= w3(y) · · ·

. . .

hr(y) =

qr−1(y) qr(y)

= wr(y)

The functions, w1(y), w2(y), · · · , wr(y), are determined from

w1(y) = h1(y) h2(y), w2(y) = h2(y) h3(y), · · · , wr−1(y) = hr−1(y) hr(y)

until

wr(y) = hr(y)

17

✬ ✫ ✩ ✪

The equations

w1(y) = 0, w2(y) = 0, · · · , wr(y) = 0

contain only simple roots, and they yield the simple, double, triple, etc. roots of f(y)

• If y0 is a root of wi(y), then it is a root of multiplicity i of f(y)
• This is a divide-and-conquer algorithm.
• The multiplicities of the roots are defined by the values of the subscripts i of the

polynomials wi(y).

• There are two essential operations

– GCD computations – Polynomial divisions but these operations are ill-posed.

18

✬ ✫ ✩ ✪ 4.1. GCD Computations

Several methods can be used to calculate the GCD of two polynomials:

• Euclid’s algorithm, The Sylvester resultant matrix, approximate polynomial fac-

torisation These methods work well if the polynomials are exact and all computations are per- formed in infinite precision arithmetic. Distinguish between the polynomials ( ˆ

f(y), ˆ g(y)) and (f(y), g(y)):

The exact polynomials

ˆ f(y) and ˆ g(y)

The noisy polynomials f(y) = ˆ

f(y) + δf(y) and g(y) = ˆ g(y) + δg(y) deg f(y) = deg ˆ f(y) = m and deg g(y) = deg ˆ g(y) = n

19

✬ ✫ ✩ ✪

• Consider an approximate greatest common divisor (AGCD) d(y) of f(y) and

g(y), f(y) ≈ u(y)d(y), g(y) ≈ v(y)d(y)

where u(y) and v(y) are coprime polynomials. Use the Sylvester resultant matrix S( ˆ

f, ˆ g) ∈ R(m+n)×(m+n) and its subresultant

matrices Sk( ˆ

f, ˆ g) ∈ R(m+n−k+1)×(m+n−2k+2) of ˆ f(y) and ˆ g(y), to calculate

the degree and coefficients of t(y)

• t = deg GCD( ˆ

f, ˆ g) = m + n − rank S( ˆ f, ˆ g)

• t = maxk rank Sk( ˆ

f, ˆ g) < m + n − 2k + 2

• The coefficients of the GCD of ˆ

f(y) and ˆ g(y) are calculated from the LU or QR

decompositions of S( ˆ

f, ˆ g).

20

✬ ✫ ✩ ✪

The coefficients of the AGCD Use the method of structured non-linear total least norm applied to the Sylvester subresultant matrix St(f, g) and the approximate polynomial factorisation of f(y) and g(y).

• Sylvester subresultant matrix St(f, g):

Calculate the structured matrix E(z) such that rank

• St(f, g) + E(z)
• = m + n − t
• Approximate polynomial factorisation of f(y) and g(y):

Calculate the polynomials p(y), q(y), r(y) and s(y) such that

f(y) + p(y) =

• u(y) + r(y)
• d(y)

g(y) + q(y) =

• v(y) + s(y)
• d(y)

21

✬ ✫ ✩ ✪ 4.2. Polynomial division

The polynomial root solver requires the computation of the polynomial h(y)

h(y) = f(y) g(y)

In practical applications, the given polynomials are inexact, and thus the computed division is

˜ h(y) = ˜ f(y) ˜ g(y) = f(y) + δf(y) g(y) + δg(y)

which yields a rational function, not a polynomial.

• How can the problem of polynomial division be posed such that the result is a

polynomial and not a rational function?

22

✬ ✫ ✩ ✪

The polynomial h(y) = f(y)/g(y), where deg f(y) = m, deg g(y) = n and

deg h(y) = m − n, is the solution of the linear algebraic equation T(g)h = f,

h =

• T(g)T T(g)

−1 T(g)T f,

where T(g) ∈ R(m+1)×(m−n+1), h ∈ Rm−n+1 and f ∈ Rm+1.

• If

r = T(g)h − f = 0 h(y) is a polynomial and g(y) is a divisor of f(y).

• If

r = T(g)h − f > 0 h(y) is a polynomial approximation of a rational function, and g(y) is not a

divisor of f(y).

23

✬ ✫ ✩ ✪

It is assumed f(y) and g(y) are corrupted by noise, and thus r > 0.

• If it is required that the result of the division f(y)/g(y) be a polynomial, repose

the problem of polynomials division as the computation of the polynomials s(y) and t(y), of minimum magnitude, such that

r(y) = f(y) + s(y) g(y) + t(y)

is a polynomial. Express this division as a matrix-vector product, such that it is necessary to com- pute the vector r such that

(T(g) + E(t)) r = f + s

where E(t) is a Tœplitz matrix whose entries are the coefficients of the polyno- mial t(y), and s is a vector whose entries are the coefficients of the polynomial

s(y).

24

✬ ✫ ✩ ✪

The equation

(T(g) + E(t)) r = f + s

has an infinite number of solutions (t, r, s).

• Since the perturbations δf(y) and δg(y) are assumed to be small, the polyno-

mials s(y) and t(y) of minimum magnitude are chosen:

min

• s2 + t2

such that

(T(g) + E(t)) r = f + s

• This is a least squares minimisation, subject to an equality constraint, the LSE
• problem. It can be solved by the QR decomposition.

25

✬ ✫ ✩ ✪

• 5. EXAMPLES

Example 5.1 εc = 10−8 exact exact computed computed relative root mult. root mult. error

• 7.5947e+00

6

• 7.59470221e+00

6 2.91588518e−07 6.3371e−01 5 6.33710738e−01 5 1.16441102e−06 1.4923e+00 5 1.49229735e+00 5 1.77821917e−06 5.4862e+00 4 5.48619195e+00 4 1.46791740e−06

• 3.3076e+00

3

• 3.30759766e+00

3 7.08230948e−07

• 3.0670e+00

2

• 3.06700222e+00

2 7.24370650e−07 4.2244e−01 2 4.22439276e−01 2 1.71431830e−06 2.5090e+00 2 2.50901393e+00 2 5.55099014e−06

26

✬ ✫ ✩ ✪

Example 5.2 εc = 10−8 exact exact computed computed relative root mult. root mult. error

• 5.9260e−01

11

• 5.92599999e−01

11 9.82710931e−10 2.6760e+00 8 2.67600001e+00 8 2.30992471e−09 6.2900e−01 5 6.28999997e−01 5 4.48781129e−09

• 9.7181e+00

4

• 9.71810000e+00

4 4.11769097e−10

27

✬ ✫ ✩ ✪

• 6. SUMMARY
• It has been shown that multiple roots of a polynomial are unstable and this makes

their accurate computation difficult.

• The pejorative manifolds of a polynomial provide a geometric interpretation of the

polynomial root finder developed by Gauss.

• The computational implementation of this root finder is not trivial because it con-

tains operations that are ill-posed.

• The method of structure non-linear total least norm is used for the solution of

these ill-posed problems.

• Computational results were presented and it was shown that structured matrix

methods are effective for the computation of multiple roots of a polynomial.

28