On a root-finding approach to the polynomial eigenvalue problem - - PowerPoint PPT Presentation

On a root-finding approach to the polynomial eigenvalue problem

Dario A. Bini

Dipartimento di Matematica, Universit` a di Pisa www.dm.unipi.it/bini

Joint work with Vanni Noferini and Meisam Sharify

Limoges, May, 10, 2012

Dario A. Bini Polynomial eigenvalue problems

The problem

Given m × m matrices Ai, i = 0, . . . , n compute the eigenvalues λ such that A(λ)v = 0, v = 0, where A(λ) =

n

• i=0

λiAi, det A(λ) ≡ 0 Standard matrix approach: linearization (A − λB)w = 0, for suitable mn × mn matrices A, B Polynomial approach: solving the scalar equation a(λ) = 0, where a(λ) = det(A(λ))

Dario A. Bini Polynomial eigenvalue problems

Structured matrix polynomials

◮ Skew-symmetric coefficients of even size: eigenvalues with

even multiplicities since a(x) = q(x)2 roots appear in pairs (λ, λ)

◮ Palindromic polynomials: Ai = An−i

T-palindromic polynomials: Ai = AT

n−i

roots appear in pairs (λ, 1/λ) including eigenvalues at infinity

◮ Hamiltonian polynomials: A2i is Hamiltonian, A2i+1 is

skew-Hamiltonian roots appear in pairs (λ, −λ)

◮ General property of the roots: roots appear in pairs (λ, f (λ)),

where f (z) is such that f (f (z)) = z Problem: to design a polynomial eigensolver which exploits the structure, and delivers approximations in pairs (λ, f (λ))

Dario A. Bini Polynomial eigenvalue problems

The Ehrlich-Aberth iteration

Given a scalar polynomial a(x) of degree N, given approximations x(0)

1 , . . . , x(0) N

to the roots λ1, . . . , λN, the Ehrlich-Aberth method generates the sequence x(ν+1)

i

= x(ν)

i

− N(x(ν)

i

) 1 − N(x(ν)

i

) N

j=1, j=i 1 x(ν)

i

−x(ν)

j

, i = 1, . . . , N where N(x) = a(x)/a′(x) is the Newton correction This iteration is just Newton’s iteration applied to the rational functions (implicit deflation) gi(x) = a(x) n

j=1, j=i(x − xj)

Dario A. Bini Polynomial eigenvalue problems

Some features of the E-A iteration

◮ Local convergence to simple roots is cubic ◮ Local convergence to simple roots is more than cubic in the

Gauss-Seidel fashion

◮ Local convergence to multiple roots is linear ◮ Global convergence is not proved: practically, convergence

always holds if initial approximations are equally placed along a circle

◮ Cost per iteration: O(N2) + N ∗ Cost of computing N(x) ◮ Number of iterations for arriving at numerical convergence is

practically bounded under a suitable choice of starting approximations

Dario A. Bini Polynomial eigenvalue problems

Problems related to the E-A implementation

◮ computing N(x) ◮ choosing initial approximations ◮ design an effective stop condition ◮ design a posteriori error bounds ◮ design effective strategies for accelerating convergence in case

• f clusters or multiple roots

These issues were treated in the package MPSolve http://www.dm.unipi.it/cluster-pages/mpsolve for approximating polynomial zeros up to any number of (guaranteed) digits [Bini, Fiorentino, Numer. Algo. 2000] We try to do the same for (structured) matrix polynomials

Dario A. Bini Polynomial eigenvalue problems

Exploiting the structure: Newton’s correction

For a general m × m matrix polynomial A(x) of degree n, set a(x) = det A(x). Then N(x) = a(x) a′(x) = 1 trace(A−1(x)A′(x)) Computational cost: O(m3 + m2n) ops

Dario A. Bini Polynomial eigenvalue problems

Exploiting the structure: Newton’s correction

Structured case (palindromic and T-palindromic polynomials): Assume that N = mn is even, set z = z(x) = x + 1/x then q(z) = x(z)−mn/2a(x(z)) is a polynomial of degree mn/2, where x(z) = (z − √ z2 − 4)/2 or x(z) = (z + √ z2 − 4)/2 is any

• f the two branches of the inverse of the function z(x).

Dario A. Bini Polynomial eigenvalue problems

Exploiting the structure: Newton’s correction

Structured case (palindromic and T-palindromic polynomials): Assume that N = mn is even, set z = z(x) = x + 1/x then q(z) = x(z)−mn/2a(x(z)) is a polynomial of degree mn/2, where x(z) = (z − √ z2 − 4)/2 or x(z) = (z + √ z2 − 4)/2 is any

• f the two branches of the inverse of the function z(x). Moreover,

q(z) q′(z) = 1 − 1/x2 g(x) − mn/(2x), g(x) = trace(A(x)−1A′(x)) This way, the computation of the roots of the palyndromic polynomial a(x) of degree mn is reduced to computing the roots of the polynomial q(z) of degree mn/2

Dario A. Bini Polynomial eigenvalue problems

Exploiting the structure: Newton’s correction

Structured case (palindromic and T-palindromic polynomials): Assume mn odd In this case −1 is root. Moreover the matrix polynomial

• A(x) = (x + 1)A(x)

is such that (n + 1)m is even and −1 has multiplicity at least m + 1. Aberth iteration can be applied to A(x) with m + 1 components of the initial vector equal to −1.

Dario A. Bini Polynomial eigenvalue problems

Exploiting the structure: Newton’s correction

Structured case: Hamiltonian polynomials (mn is even) Set z = x2. For x(z) = √z or x(z) = −√z, q(z) = a(x(z)) is a polynomial of degree mn/2. Moreover, q(z) q′(z) = 2x g(x) − 1

x

, g(x) = trace(A(x)−1A′(x)) If mn is odd, then q(z) = a(x(z))/x(z) is a polynomial and q(z) q′(z) = 2x3 g(x) − 1, g(x) = trace(A(x)−1A′(x)) Once the roots z1, . . . , zmn/2 have been computed, the roots of a(x) are given by the pairs (√zi, −√zi), i = 1, . . . , mn/2

Dario A. Bini Polynomial eigenvalue problems

Exploiting the structure: Newton’s correction

General case: roots in pairs (x, f (x)), f (x) = (αx + β)/(γx + α), α2 + βγ = 0 Set z(x) = xf (x) = αx2+βx

γx−α

and denote x(z) one of the two branches of the inverse of z(x). If there are no eigenvalues with

• dd multiplicity, then

q(z) = a(x(z)) (γx(z) − α)mn/2 is a polynomial

Dario A. Bini Polynomial eigenvalue problems

Computational cost

The overall cost of E-A iteration applied to an m × m matrix polynomial of degree n is O(m4n + m3n2) ops The overall cost of the technique based on linearization is O(m3n3) ops From the complexity point of view the E-A approach is more convenient if n > m The cost of E-A, i.e., the computation of trace(A(x)−1A′(x)), can be reduced to O(m2n2) ops, by using linearization. This leads to the overall cost of O(m3n3) ops The drawback is that linearization may increase the condition number of eigenvalues.

Dario A. Bini Polynomial eigenvalue problems

Choosing initial approximations

The choice of initial approximations is crucial to arrive at numerical convergence with a small number of iterations This is particularly important for polynomials having very unbalanced roots The package MPSolve for computing polynomial roots to any guaranteed precision, relies on a very effective tool: the Newton polygon construction We synthesize the main ideas of this tool and then extend it to the case of matrix polynomials.

Dario A. Bini Polynomial eigenvalue problems

Choosing initial approximations

The Rouch´ e theorem: If p(x) and q(x) are polynomials such that |p(x)| > |q(x)| for |x| = r then p(x) and p(x) + q(x) have the same number of roots in the

• pen disk of center 0 and radius r

Given a(x) = n

i=0 aixi apply Rouch´

e theorem with p(x) = akxk and q(x) = n

i=0, i=k xiai. It follows that if

|p(x)| = |ak||x|k >

n

• i=0,i=k

|ai||x|i ≥ |q(x)|, |x| = r then a(x) = p(x) + q(x) has k roots of modulus less than r.

Dario A. Bini Polynomial eigenvalue problems

Choosing initial approximations

Properties: The equation |ak|rk =

n

• i=0,i=k

|ai|ri has

◮ one positive solution t0 if k = 0 ◮ one positive solution sn if k = n ◮ either no positive solutions or two positive solutions s ≤ t

• therwise

Pellet’s Theorem, 1881 If h1, . . . , hp are the values of k for which the above equation has either one or two positive solutions shi ≤ thi then

◮ the open annulus with radii shi, thi contains no roots of any

polynomial equimodular with p(x)

◮ the closed annulus with radii thi−1, shi contains hi − hi−1 roots

Dario A. Bini Polynomial eigenvalue problems

Choosing initial approximations

Dario A. Bini Polynomial eigenvalue problems

Choosing initial approximations: strategy of MPSolve

Computing the positive roots of n polynomials is expensive. An alternative solution relies on the following properties:

◮ Any positive solution r of the blue equation, if it exists, is

such that uk := max

i<k

• ai

ak

• 1

k−i

< r < min

i>k

• ai

ak

• 1

k−i

=: vk

◮ Let ki, i = 1, . . . q be such that uki < vki. Then vki = uki+1,

moreover any interval [shj, thj] does not contain any uki

◮ The indices ki are the abscissas of the Newton polygon: upper

convex hull of the set {(i, log |ai|), i = 0, . . . , n}

◮ The values uki are the exponential of the slopes of the edges

• f the Newton polygon. They are called tropical roots

Dario A. Bini Polynomial eigenvalue problems

Choosing initial approximations: strategy of MPSolve

Dario A. Bini Polynomial eigenvalue problems

Choosing initial approximations: strategy of MPSolve

Dario A. Bini Polynomial eigenvalue problems

Starting approximations: an example

p(x) = x10 + (1000x + 1)3 Coefficient vector: (1, 3000, 3000000, 1000000000, 0, 0, 0, 0, 0, 0, 1) There are: 3 roots of modulus close to 1/1000, 7 roots of modulus close to 19.3065 The Newton polygon has three vertices besides 0 and n. The values of uk are 0.0003, 0.001, 0.003, 19.3069 7 initial approximations are placed in a circle of radius 19.3069 3 initial approximations are placed inside a disk of radius 3/1000

Dario A. Bini Polynomial eigenvalue problems

The Sharify improvement

What happens inside a blue annulus? If a vertex of the polynomial is sufficiently sharp then we may arrive at a strong localization of the roots Denote

◮ ri = uki the value of the radius corresponding to the vertex ki ◮ mi = ki+1 − ki its “multiplicity” ◮ δi = ri/ri+1 the ratio of two consecutive radii

Theorem (Sharify) If δi, δi−1 < 1/9 then the polynomial a(x) has mi roots in the annulus with radii ri/3 and 3ri.

Dario A. Bini Polynomial eigenvalue problems

Extension to matrix polynomials

Let A, B Hermitian matrices, define A ≻ B if A − B is positive definite, and A B if A − B is positive semidefinite Let F(x), G(x) be matrix polynomials

Theorem (Rouch´ e theorem for matrix polynomials)

If F(x)∗F(x) ≻ G(x)∗G(x) for |x| = r then F(x) and F(x) + G(x) have the same number of roots in the disc of center 0 and radius r.

Corollary

If (A∗

kAk)r2k ≻ (n i=0, i=k Aixi)∗(n i=0, i=k Aixi),

for |x| = r then det A(x) has kn roots inside the disk of center 0 and radius r. The above condition requires that det Ak = 0 and is implied by rk ≻

n

• i=0, i=k

A−1

k Airi

Dario A. Bini Polynomial eigenvalue problems

Generalization

Theorem

• 1. If 0 < k < n is such that det Ak = 0 then the equation

r k =

• i=0, i=k

A−1

k Air i

(1) has either no real positive solution or two real positive solutions sk ≤ tk.

• 2. In the latter case, the polynomial a(x) = det A(x) has no roots in

the open annulus of radii sk, tk, while it has mk roots in the disk of radius sk.

• 3. If k = 0 and det A0 = 0 then (1) has only one real positive root t0,

moreover, the polynomial a(x) has no root in the open disk of radius t0.

• 4. If k = n and det An = 0, then (1) has only one real positive solution

sn and the polynomial a(x) has no root outside the disk of radius sn.

Dario A. Bini Polynomial eigenvalue problems

Corollary (Generalized Pellet Thm.)

Let h0 < h1 < ... < hq be the values of k such that det Ak = 0 and there exist positive real solution(s) shi ≤ thi of (1). Then

• 1. thi−1 ≤ shi, i = 1, . . . , q;
• 2. there are m(hi − hi−1) roots of a(x) in the annulus of radii

thi−1, shi;

• 3. there are no roots of the polynomial a(x) in the open annulus
• f radii shi, thi, where i = 0, 1, . . . , q and we assume that

s0 = 0, tn = ∞.

Dario A. Bini Polynomial eigenvalue problems

Theorem

If k is such that (1) has two real positive solutions sk ≤ tk, then uk ≤ sk ≤ tk ≤ vk where uk := maxi<k ||A−1

k Ai||1/(k−i),

vk := mini>k ||A−1

k Ai||1/(k−i)

If, for k = 0 (1) has a solution t0 then t0 ≤ v0 := min

i>0 ||A−1 0 Ai||−1/i

If for k = n (1) has a solution sn then sn ≥ un := max

i<n ||A−1 n Ai||1/(n−i)

Moreover, vki ≤ uki+1, i = 1, . . . , p − 1.

Dario A. Bini Polynomial eigenvalue problems

A special class

Consider the class Qm,n of matrix polynomials A(x) = n

i=0 xiAi,

where Ai = αiQi, Q∗

i Qi = I,

αi ≥ 0 The condition given in the blue equation (1) turns into αkrk =

n

• i=0,i=k

αiri and leads to the machinery valid in the scalar case for the location

• f the roots of a polynomial by means of the Newton polygon.

This is the best result that we can obtain for polynomials in the class Qm,n. In fact, for Qi = I we obtain exactly n copies of the scalar polynomial p(x) = n

i=0 αixi for which the inclusion result

is optimal.

Dario A. Bini Polynomial eigenvalue problems

Theorem

Let S = {hi, i = 0, ..., p} be such that uk < vk if and only if k ∈ S. Then,

• 1. {k0, ..., kq} ⊂ {h0, ..., hp}, moreover

{h0, ..., hp} ∩ [ski, tki] = ∅;

• 2. vki ≤ uki+1;
• 3. if A(x) ∈ Qm,n, then vki = uki+1 and vki coincide with the

vertices of the Newton polygon of the polynomial w(x) = n

i=0 Aixi = n i=0 αixi.

Dario A. Bini Polynomial eigenvalue problems

The Sharify bound for matrix polynomials

Let A(x) = n

i=0 xiAi ∈ Qm,n, i.e., Ai ∈ Cm×m, Ai = αiQi,

Q∗

i Qi = I. Let a(x) = det A(x).

Define ri the radii (tropical roots) corresponding to the i-th vertex

• f the Newton polygon of the polynomial

n

• i=0

xiAi2 and mi their multiplicities, define δi = ri/ri−1 Theorem For A(x) ∈ Qm,n, if δi, δi−1 < 1/12.12 then the polynomial a(x) det A(x) has mmi roots in the annulus with radii ri/4.37 and 4.37ri.

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Counting # of iterations

Orthogonal random coefficients, degree n = 13, scaling factors: α = [1, 3 · 103, 3 · 106, 109, 0, . . . 0, 104, 0, 0, 1] Newton polygon Unit circle m simul it aver it simul it aver it 5 8 5.4 243 191 10 9 5.5 444 375 20 11 5.6 855 738 40 13 6.1 1594 1466

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Counting # of iterations

Non-orthogonal random coefficients, degree n = 13, scaling factors: α = [1, 3 · 103, 3 · 106, 109, 0, . . . 0, 104, 0, 0, 1] Newton polygon Unit circle m simul it aver it simul it aver it 5 9 6.8 240 190 10 13 7.7 457 372 20 16 9 851 732 40 16 10.4 1597 1457

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Counting # of iterations

General random coefficients, scaling factors: coefficient i, factor: (rand < 1/4) ∗ 1012 Number of average iterations with m = 10 and n = 30 Newton’s Polygon 7.2 Unit circle 111

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Approximation error

Ehrlich-Aberth: ∗ Polyeig: + Problem: Sommerfeld problem

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Approximation error

Ehrlich-Aberth: ∗ Polyeig: + Problem: Plasma drift

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Approximation error

Ehrlich-Aberth: ∗ Polyeig: + QuadEig (Hammarling-Munro-Tisseur): × Problem: hospital

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Approximation error

Ehrlich-Aberth: ∗ Polyeig: + QuadEig (Hammarling-Munro-Tisseur): × Problem: power plant

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Approximation error

Structured Ehrlich-Aberth: Ehrlich-Aberth: ∗ Polyeig (QZ): + Structured URV (Schroeder): × Problem: butterfly

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Approximation error

Structured Ehrlich-Aberth: Ehrlich-Aberth: ∗ Polyeig (QZ): + Structured URV (Schroeder): × Problem: wiresaw1

Dario A. Bini Polynomial eigenvalue problems

Some numerical experiments: Approximation error

Structured Ehrlich-Aberth: Ehrlich-Aberth: ∗ Polyeig (QZ): + Problem with symmetry (λ, λ+1

λ−1)

Dario A. Bini Polynomial eigenvalue problems

Some references

–O. Aberth. Iteration methods for finding all zeros of a polynomial simultaneously.

• Math. Comput, 1973.

–D. A. Bini and G. Fiorentino. Design, analysis, and implementation of a multiprecision polynomial rootfinder. Numer. Algorithms, 2000. –D.A. Bini, V. Noferini. Solving polynomial eigenvalue problems by means of the Ehrlich-Aberth method. TR 2012. –L. Gemignani and V. Noferini. The Ehrlich-Aberth method for palindromic matrix polynomials represented in the Dickson basis. Linear Algebra Appl. 2011 –H. Guggenheimer, Initial Approximations in Durand-Kerner’s Root Finding Method, BIT, 26 (1986). –S. Hammarling, C. Munro and F. Tisseur. An Algorithm for the Com- plete Solution

• f Quadratic Eigenvalue Problems. MIMS EPrint 2011.86, avaliable online.

–N. J. Higham, D. S. Mackey and F. Tisseur. The conditioning of linearizations of matrix polynomials. SIAM J. Matrix Anal. Appl., 2006. –V. Noferini. The application of the Ehrlich-Aberth method to structured polynomial eigenvalue problems. Proc. Appl. Math. Mech., 2011. –V. Noferini. Polynomial Eigenproblems: A Root-Finding Approach, Ph.D. Thesis, 2012, Dipartimento di Matematica Universit` a di Pisa, Italy. –A. Ostrowski, On a Theorem by J.L. Walsh Concerning the Moduli of Roots of Algebraic Equations, Bull. A.M.S., 47 (1941) 742–746.

Dario A. Bini Polynomial eigenvalue problems

Some references

–A.E. Pellet, Sur un mode de s´ eparation des racines des ´ equations et la formule de Lagrange, Bulletin des Sciences Math´ ematiques, (2), vol 5 (1881), pp. 393–395. –M. Sharify. Scaling Algorithms and Tropical Methods in Numerical Ma- trix Analysis: Application to the Optimal Assignment Problem and to the Accurate Computation of

• Eigenvalues. Ph.D. Thesis, Ec´
• le Polytechnique Paris Tech, 2011

–C. Schr¨

• der. Palindromic and Even Eigenvalue Problems – Analysis and Numerical
• Methods. PhD thesis, T.U. Berlin, 2008.

– P.P. Vaidyanathan, S.K. Mitra, A unified structural interpretation of some well-known stability test procedures for linear systems. Proceedings of the IEEE, vol 75, no. 4, 1987. –J.L. Walsh, On Pellet’s theorem concerning the roots of a polynomial, Annals of Mathematics, vol 26 (1924), pp. 59–64.

Dario A. Bini Polynomial eigenvalue problems