Tropical roots as approximations to eigenvalues of matrix polynomials

The tropical roots of t×p(x) = max0≤j≤` ‖Aj‖x are points at which the maximum is attained at least twice. These roots, which can be computed in only O(`) operations, can be good approximations to the moduli of the eigenvalues of the matrix polynomial P (λ) = ∑` j=0 λ Aj , in particular when the norms of the matrices Aj vary widely. Our aim is to investigate this observation and its applications. We start by providing annuli defined in terms of the tropical roots of t×p(x) that contain the eigenvalues of P (λ). Our localization results yield conditions under which tropical roots offer order of magnitude approximations to the moduli of the eigenvalues of P (λ). Our tropical localization of eigenvalues are less tight than eigenvalue localization results derived from a generalized matrix version of Pellet’s theorem but they are easier to interpret. Tropical roots are already used to determine the starting points for matrix polynomial eigensolvers based on scalar polynomial root solvers such as the Ehrlich-Aberth method and our results further justify this choice. Our results provide the basis for analyzing the effect of Gaubert and Sharify’s tropical scalings for P (λ) on (a) the conditioning of linearizations of tropically scaled P (λ) and (b) the backward stability of eigensolvers based on linearizations of tropically scaled P (λ). We anticipate that the tropical roots of t×p(x), on which the tropical scalings are based, will help designing polynomial eigensolvers with better numerical properties than standard algorithms for polynomial eigenvalue problems such as that implemented in the MATLAB function polyeig.


Introduction.
Being able to cheaply locate the eigenvalues of a real or complex n × n matrix polynomial (1.1) is useful in a number of situations, such as, for example, when selecting the starting points in the Ehrlich-Aberth method for the numerical solution of polynomial eigenvalue problems [6], [7], or in choosing the contour in contour integral methods for polynomial eigenvalue problems of large dimensions [3].Betcke's diagonal scaling [4, section 5], whose aim is to improve the conditioning of P 's eigenvalues near a target eigenvalue ω, requires a priori knowledge of the magnitude of ω.
The tropical roots of a tropical (or max-times) polynomial f (x) = max 0≤i≤ (a i x i ) with a i , x ≥ 0 are points (i.e., nonnegative real numbers) at which the maximum is attained for at least two values of i for some x.They are easy and cheap to compute (see section 2.1).Our aim is to investigate the order of magnitude approximation of the eigenvalues of P (λ) in terms of the tropical roots of t × p(x) = max 0≤i≤ ( A i x i ) for some matrix norm • subordinate to a vector norm.
Gaubert and Sharify [9,Thm. 2] were the first to notice the tropical splitting of the eigenvalues of matrix polynomials.Indeed, for n × n heavily damped quadratics, i.e., quadratic matrix polynomials Q(λ) = λ 2 A 2 + λA 1 + A 0 with A 1 2 ≥ A 0 A 2 , they showed that gap Λ(Q), Λ(L) ≤ g(κ(A 2 ))α max α min α max 1/(2n) , (1.2) (1.3) where Λ(P ) denotes the spectrum of P (λ), gap(Λ(Q), Λ(L)) is a measure of the distance between the n largest eigenvalues of Q(λ) in modulus, and the n eigenvalues of L(λ) = A 2 λ + A 1 , g(κ(A 2 )) is more or less a constant times the matrix condition number κ(A 2 ) = A 2 A −1 2 , and α max and α min are the largest and smallest tropical roots of t × q(x) := max( A 0 , A 1 x, A 2 x 2 ).The bounds (1.2)- (1.3) show that when the ratio α min /α max is small enough and A 2 , A 1 are well conditioned then there are precisely n eigenvalues of Q(λ) with moduli of the order of α max .Similarly, when A 1 and A 0 are both well conditioned, the moduli of the n smallest eigenvalues of Q(λ) are close to the smallest tropical root α min of t × q(x).
For the particular case of matrix polynomials P (λ) with coefficient matrices of the form A i = σ i Q i with σ i ≥ 0 and Q * i Q i = I, and the 2-norm • 2 , Bini, Noferini, and Sharify [7,Thm. 2.7] have identified annuli of small width defined in terms of the tropical roots of t × p(x) that contain the eigenvalues of P (λ).We extend their results to arbitrary matrix polynomials and any subordinate matrix norm in section 2.2.We obtain conditions under which tropical roots offer order of magnitude approximations to the moduli of the eigenvalues of P (λ).As shown in section 3 our tropical localization results are less tight than those from the generalized Pellet's theorem, both in the form given in [7,Thm. 2.1] and in [17,Thm. 3.3], but they are easier to interpret.We illustrate our localization results with numerical examples in section 4 and show experimentally how tropical roots can help in the design of a numerically stable polynomial eigensolver.
We note that a different approach, also involving tropical roots, is pursued in [2], where Akian, Gaubert, and Sharify derive bounds for products of eigenvalues of n × n matrix polynomials P (λ) of degree .Their results generalize to matrix polynomials bounds by Ostrowski [18,19] for products of roots of scalar polynomials.

Tropical bounds.
The max-plus semiring R max is the set R∪{−∞} equipped with the max operation denoted by ⊕ as addition and the usual addition denoted by ⊗ as multiplication.The zero and unit elements of this semiring are −∞ and 0, respectively.
A variant of R max is the max-times semiring R max,× , which is the set of nonnegative real numbers R + equipped with the max operation as addition and the usual multiplication as multiplication.This semiring is isomorphic to R max by the map x → log x.So, every notion defined over R max has an R max,× analogue.By the word "tropical," we refer to any of these algebraic structures.
Fig. 1.Newton polygon (in blue) corresponding to a max-plus tropical polynomial tp(x) = j=0 a j ⊗ x ⊗j .Points (j, a j ) are denoted by red dots.

Tropical polynomial and Newton polygon.
A max-plus tropical polynomial tp is a function of a variable x ∈ R max of the form where is a nonnegative integer and a 0 , . . ., a ∈ R max .The tropical polynomial tp is of degree if a = −∞.If we assume that at least one of the coefficients a 0 , . . ., a is finite then tp is a real valued convex function, piecewise affine, with integer slopes.The finite tropical roots of tp(x) are the points at which the maximum in the expression (2.1) is attained for at least two values of i for some x.The tropical roots can be obtained via Newton polygons.Define the Newton polygon of tp to be the upper boundary of the convex hull of the set of points (j, a j ), j = 0, . . ., (see Figure 1).This boundary consists of a number of linear segments.The opposites of the slopes of these segments are precisely the tropical roots and the multiplicity of a root coincides with the width of the corresponding segment measured by the difference of the abscissae of its endpoints (see [1,Prop. 2.10] or [16,Lem. 2.3]).Hence, if we denote by k 0 = 0 < • • • < k q = the abscissae of the vertices of the Newton polygon then tp(x) has q distinct roots given by (2.2) with multiplicities m j = k j − k j−1 , j = 1, . . ., q, respectively.Since the points (j, a j ) are already sorted by abscissae, the Graham scan algorithm [12] computes the convex hull of these points in O( ) operations.As a result the tropical roots, counted with multiplicities, can be computed in O( ) operations [9, Prop.1].
In the "max-times" semiring R max,× , a tropical polynomial has the form t × p(x) = max 0≤i≤ a i x i , where a 0 , . . ., a are nonnegative numbers, and x takes nonnegative values.The tropical roots of t × p(x) are, by definition, the exponentials of the tropical roots of the max-plus polynomial tp(x) = max 0≤i≤ (log a i + ix).In view of (2.2), the q distinct tropical roots of t × p(x) and their multiplicities are given by (a kj−1 /a kj ) 1/(kj −kj −1) , m j = k j − k j−1 , j = 1, . . ., q, respectively.

Eigenvalue location: Tropical approach.
Throughout the rest of this paper, any matrix polynomial denoted by P (λ) is regular, i.e., det P (λ) is not identically zero.Moreover, we assume for simplicity that A 0 = 0. Note that this is no loss of generality for the purpose of eigenvalue location, since otherwise, if we let κ := min{i s.t.A i = 0}, then there must be at least nκ zero eigenvalues and we may simply consider λ −κ P (λ).Then the finite eigenvalues of an n × n P (λ) of degree are the roots of det P (λ) = 0, and if det P (λ) = 0 has degree d ≤ n, then P (λ) has n − d eigenvalues at infinity.To P (λ) = i=0 A i λ i we associate the max-times tropical scalar polynomial where • is any matrix norm subordinate to a vector norm.Our aim is to show that, under specific conditions, the tropical roots of t × p(x) are good approximations to the moduli of the eigenvalues of P (λ).The key tool for this is a generalization of Rouché's theorem for matrix valued functions [10], [17].Theorem 2.1 (generalized Rouché theorem for matrix valued functions).Let P, Q : Ω → C n×n be analytic matrix valued functions, where Ω is an open connected subset of C and assume that P (λ) is nonsingular for all λ on the simple closed curve Γ ⊆ Ω.Let • be any matrix norm on C n×n induced by a vector norm on C n .If P (λ) −1 Q(λ) < 1 for all λ ∈ Γ, then det(P + Q) and det(P ) have the same number of zeros inside Γ, counting multiplicities.
Before deriving our new results, we set up the notation used throughout the reminder of this paper.

Notation.
The variable i will usually be an index varying between 0 and the degree of t × p(x) in (2.3), whereas j will be an index with value between 1 and q, where q is the number of distinct tropical roots of t × p(x).These tropical roots will be denoted by α j , j = 1, . . ., q, with (2.4) where denote the abscissae of the Newton polygon associated with the max-plus polynomial tp(x) = max 0≤i≤ (log A i + ix) (see Figure 1).We write (2.5) Note that the tropical roots (2.4) have the property that 0 As in [7] and [9], the tropical roots (2.4) will be used to define an eigenvalue parameter scaling, λ = α j μ, and a scaled matrix polynomial, P (μ), via (2.6) t × p(α j ) Clearly, μ is an eigenvalue of P (μ) if and only if α j μ is an eigenvalue of P (λ).Note that this scaling does not affect the condition number with respect to inversion, We will use disks and annuli to localize the eigenvalues of P (λ).The closed (open) disk in the complex plane centered at 0 with radius r is denoted by D(r) ( the ratio between two consecutive roots.

Preliminary results.
We now present preliminary lemmas, which will be needed in section 2.2.3 to prove our localization results.We refer to section 2.2.1 for the notation.
The norms of the coefficient matrices of the scaled matrix polynomial P (μ) in (2.6) are at most 1 as shown by this first lemma.
The next lemma provides upper and lower bounds on the moduli of all the eigenvalues of P (λ) in terms of the smallest and largest tropical roots α 1 and α q of t × p(x), and the conditioning of A 0 and A .
Lemma 2.3.Let P (λ) = i=0 A i λ i with A 0 , A = 0 be a regular matrix polynomial.Then every eigenvalue of P (λ) satisfies Furthermore, if both A 0 and A are invertible, both inequalities are strict.Proof.For the upper bound, we consider the scaled matrix polynomial P (μ) in (2.6) with j = q.Observe first that if A is singular then the right-hand side is ∞, so there is nothing to prove and, if P (λ) is regular, then the bound is attained in the sense that necessarily P (λ) has an eigenvalue at infinity.Hence, we may assume that A is invertible.Let θ = max 0≤i≤ −1 A i 1/( −i) .We now recall an argument from Hence, any eigenvalue must satisfy |μ| ≤ β := (1 + A −1 )θ (observe that θ < β so this does not contradict the assumption |μ| > θ).By Lemma 2.2, and since A 0 = 0, we have θ = 1 > 0, while κ(A ) < ∞ implies |μ| < ∞.Thus, the argument above can be tightened as the inequality i=1 The lower bound is proved similarly, using P (μ) in (2.6) with j = 1 and applying the above argument to rev P (μ) = μ P (1/μ).We invite the reader to fill in the details.
With the aim of invoking Theorem 2.1, we decompose P (μ) = (t × p(α j )) −1 P (λ) as the sum of two matrix polynomials, (2.9) We will need the following localization result for the nonzero eigenvalues of S(μ).Lemma 2.4.If A kj−1 and A kj are nonsingular then the nm j nonzero eigenvalues of S(μ) in (2.9) are located in the open annulus Proof.Note that S(μ) is regular, of degree k j , with nk j−1 zero eigenvalues.Hence S(μ) has nk j −nk j−1 = nm j nonzero eigenvalues, which are eigenvalues of μ −kj−1 S(μ).The tropical polynomial associated with μ −kj−1 S(μ) has only one root, which is equal to 1.The lemma is then a direct consequence of Lemma 2.3.
We now consider μ ∈ C such that |μ| ≥ 1 + κ(A kj ).Note that such a μ exists only if A kj is nonsingular.Lemma 2.4 implies that for such a μ, the matrix S(μ) is invertible.We rewrite S(μ) as in (2.10) with The rest of the proof is then analogous to the case where μ ∈ S so we omit it.
Finally, this last technical lemma will be needed in the proof of our tropical localization results in section 2.2.3, and in section 3 when comparing the Pellet bounds with the tropical bounds.
Lemma 2.6.For given c, δ > 0 such that δ ≤ (1+2c) −2 , the quadratic polynomial (i) The discriminant of p(r) is not negative since δ ≤ (1 + 2c) −2 so p has two real roots, f, g such that f ≤ g.That f ≥ 1 + c is easy to check.(ii) Clearly, Since f and g are the roots of p, f

Main results. When
for all 0 ≤ i ≤ then P (λ) has only one tropical root given by α = ( A 0 / A ) 1/ and we know from Lemma 2.3 that all the eigenvalues of P (λ) lie in the annulus So in this case, if A 0 and A are well conditioned then P (λ) has n eigenvalues of modulus close to α.We now extend this type of result to the case where t × p(x) has more than one tropical root.Theorem 2.7.
, where f (δ, c) is defined as in Lemma 2.6, and g j = (δ j f j ) −1 .Then, in the notation of section 2.2.1, the following statements hold.

and it does not have any eigenvalue in the open annulus
• A(f j α j , g j α j ).(ii) If δ j ≤ (1 + 2 κ(A kj )) −2 and δ s ≤ (1 + 2 κ(A ks )) −2 with 1 < j < s < q then P (λ) has exactly n(k s − k j ) eigenvalues inside the annulus A(g j α j , f s α s ).(iii) Let j be the smallest index, if any, such that δ j ≤ (1+2 κ(A kj )) −2 .Then P (λ) has exactly nk j eigenvalues inside the annulus A((1 + κ(A 0 )) −1 α 1 , f j α j ).(iv) Let j be the largest index, if any, such that δ j ≤ (1 + 2 κ(A kj )) −2 .Then P (λ) has exactly n( − k j ) eigenvalues inside the annulus A(g j α j , (1 + κ(A ))α q ).Proof.(i) We assume that δ j ≤ (1 + 2 κ(A kj )) −2 and we partition P (μ) as in (2.9).Let r be such that Note that such r exists since we can apply the bounds in Lemma 2.5.These yield The latter bound is less than 1 if or, equivalently, if , or equivalently when p(r) < 0, where p is as in Lemma 2.6 with δ = δ j and c = κ(A kj ).It follows from Lemma 2.6 that p(r) is negative for the values of r such that (2.15) f j < r < g j recalling that by Lemma 2.6 f j and g j are the two roots of p.Note that, by the same lemma, A(f j , g j ) and has exactly nk j eigenvalues inside the disk D(f j ).This completes the proof of (i) since λ = μα j .(ii) If δ s ≤ (1 + 2 κ(A ks )) −2 then by (i), P (λ) has nk s eigenvalues inside the disk D(f s α s ) and no eigenvalues inside the annulus ).An analogous statement holds if we assume that δ j < (1 + 2 κ(A kj )) −2 .This implies that P (λ) has exactly n(k s −k j ) eigenvalues which lie in the annulus A(g j α j , f s α s ).(iii) If δ j ≤ (1 + 2κ(A kj )) −2 then by (i) P (λ) has nk j eigenvalues inside the disk D(f j α j ).Also by Lemma 2.3 P (λ) does not have any eigenvalue inside the disk ). (iv) The statement follows from (i) and Lemma 2.3.Note that for a fixed value of c ≥ 1 and δ ≤ (1 + 2c) −2 , f (δ, c) in Lemma 2.6 is an increasing function of δ and its maximum value, which is 1 + 2c, is achieved at δ = (1+2c) −2 .This implies that f (δ j , κ(A kj )) ≤ 1+2 κ(A kj ) for any δ j ≤ (1+2 κ(A kj )) −2 and we have the following corollary.
Corollary 2.8.In the notation of section 2.2.1 and under the assumptions of Theorem 2.7, the following statements hold.

Comparisons with Pellet's bounds.
Pellet's theorem, generalized to matrices by Bini, Noferini, and Sharify [7, Thm.2.1] and Melman [17,Thm. 3.3], also provides annuli containing the eigenvalues of P (λ).Although [7, Thm.2.1] is stated for the 2-norm, as mentioned in [17] the result can be extended to any subordinate norm by using Theorem 2.1.Throughout this section, roots of polynomials are counted with multiplicity: in particular, a double positive root is thought of as two coincident positive roots.
Theorem 3.1 (generalized Pellet theorem).Let P (λ) = i=0 A i λ i with > 1 and A 0 , A = 0, and let • denotes any subordinate matrix norm.For A k nonsingular define The following statements hold.).(iii) If A is nonsingular then q (x) has only one real positive root s and all the eigenvalues of P (λ) are located in the disk D(s ).Remark 3.2.Theorem 3.1 can be restated in terms of annuli in an analogous way to the statement of Theorem 2.7 and Corollary 2.8.Indeed, if A 0 (resp.,A ) is nonsingular let t 0 (resp., s ) be the unique real root of q 0 (x) (resp., q (x)), or t 0 = 0 (resp., s = 0) otherwise.Then the following statements hold. ( Although the coefficients of q k (x) are less expensive to compute than those of q k (x), Theorem 3.1 provides tighter localization results than those in [17,Thm. 3.3].This fact has appeared, in the form of a remark, in [7, p. 1713], [17, p. 1555].The next proposition provides a detailed statement.Proposition 3.3.Let P (λ) = i=0 A i λ i with > 1 and A 0 = 0. (i) If A k is nonsingular for some k such that 1 ≤ k ≤ − 1, and q k (x) has exactly two real positive roots s k , t k with s k ≤ t k then q k (x) has also exactly two real positive roots, s k , t k such that s k ≤ s k ≤ t k ≤ t k .(ii) If A 0 is nonsingular and q 0 (x) has a real positive root t 0 then q 0 (x) has also a real positive root t 0 such that t 0 ≤ t 0 .(iii) If A is nonsingular and q (x) has a real positive root s n then q k (x) has also a positive root s n such that s n ≤ s n .Proof.(i) Since for any subordinate matrix norm AB ≤ A B , and q k ( s k ) = q k ( t k ) = 0, we have which implies that q k ( s k ) ≤ 0. Similarly we can show q k ( t k ) ≤ 0. Since q k (0) > 0 and the leading coefficient of q k (x) is positive, this implies that q k (x) has exactly 1 two positive roots s k ≤ t k .But q k ( s k ), q k ( t k ) ≤ 0, so we must have The statements (ii) and (iii) are proved in a similar way.The next result, when combined with Proposition 3.3, shows that the eigenvalue localization results from either form of the generalized Pellet theorem are always better than those presented in Theorem 2.7.
Proposition 3.4.Using the notation of Theorem 2.7, assume that δ j ≤ (1 + 2 κ(A kj )) −2 for some 1 ≤ j ≤ q − 1.Then the k j th Pellet polynomial q kj (x) in (3.2) has exactly two positive roots s kj < t kj such that s kj < f j α j and t kj > g j α j .
Proof.Note that q kj (0) = A 0 ≥ 0 and lim x→∞ q kj (x) ∼ x A > 0. Also, the condition on δ j in the statement implies that A kj is nonsingular.According to Pellet's theorem, q kj has either zero or two positive roots.Hence, if there exist two positive numbers y 1 < y 2 such that q kj (y 1 ) < 0 and q kj (y 2 ) < 0, then q kj has two positive roots s kj , t kj such that s kj < y 1 < y 2 < t kj .Define y 1 := f j α j and y 2 := g j α j .Next we show that q kj (y 1 ), q kj (y 2 ) < 0. Note that since by Lemma 2.2 we have that recalling that g j = (δ j f j ) −1 .This yields , which implies that q kj (y 1 ) < 0. The proof for y 2 is similar to the one given above so we skip it.
Let q k be as in (3.1) and similarly to H but with q k (x) in place of q k (x).
Proof.The first two inclusions follow directly from Propositions 3.4 and 3.3.For the last inclusion, suppose that k ∈ K. Then there are k a < k < k b such that k a and k b are two consecutive indices in K. Modulo the appropriate scaling we may assume where we used AB ≥ B A −1 , which holds for any invertible A and any subordinate matrix norm.Therefore, if 0 < x < 1, q k (x) We conclude that q k (x) > 0 for all x > 0. So q k (x) does not have a real positive root and hence k ∈ H. Theorem 3.5 shows that to compute the Pellet bounds we only need to construct the polynomials q k (x) and q k (x) for k ∈ K.
Let s 0 = 0 and let t 0 be the real positive root of q 0 (x) if A 0 is nonsingular and t 0 = 0 otherwise.Also, let t = +∞ and let s to be the real positive root of q (x) if A is nonsingular and set s = +∞ otherwise.We define s 0 , t 0 , s , and t similarly with respect to q k (x), k = 0, .It is shown in [7, Cor.2.2] that t hj ≤ s hj+1 and it follows from Proposition 3.3 that t hj ≤ s hj+1 .Then from Theorem 3.1 and Propositions 3.3-3.4and (2.18) it follows that (3.5) with p ≥ r ≥ m.In other words, (3.5) means that Bini et al.'s generalized Pellet theorem provides better eigenvalue localization results than Melman's generalized Pellet theorem, which in turn provides better localization results than our localization theorem based on tropical roots (see Theorem 2.7 and Corollary 2.8).However, the tropical roots and the results of Theorem 2.7 and its corollary remain interesting since these results can be easily interpreted and can be used in the numerical computation of the eigenvalues as we explain below.Importantly, the amount of information that Theorem 2.7 provides does not depend on the condition numbers of all the coefficients, but only on a selected number of them (A k such that k ∈ K).This fact can be used to give bounds on the sensitivity of the moduli of the eigenvalues when one coefficient A i , i ∈ K, is perturbed.Even when the conditions of Theorem 2.7 are not satisfied, it can still happen that the tropical roots provide good approximations to the moduli of the eigenvalues.Indeed, they always lie inside the inclusion annuli defined by the generalized Pellet theorem, as we now show.
Theorem 3.6.Let P (λ) = i=0 A i λ i be a regular matrix polynomial.Also, for some 0 ≤ j ≤ p and some 0 ≤ i 1 < i 2 ≤ q, let h j = k i1 and h j+1 = k i2 be two consecutive indices in H ⊆ K, defined as in (3.3) and (2.5).Then where s hj ≤ t hj are the two positive real roots of q hj (x) in (3.1).Norm and condition number of the coefficient matrices of P (λ) in Experiment 1.

Numerical experiments and applications.
We start with some experiments that illustrate the bounds of sections 2.2 and 3 and show how well the tropical roots of t × (x) approximate the moduli of the eigenvalues of P (λ).Our experiments were performed in MATLAB 7, for which the unit roundoff is u = 2 −53 ≈ 1.1 × 10 −16 .
Experiment 1.Our first example is a 20 × 20 quartic matrix polynomial P (λ) = 4 i=0 λ i A i generated with the MATLAB commands randn('state',48); n = 20; A0 = 1e-5*randn(n); A1 = 1e2*randn(n); A2 = 1e2*randn(n); A3 = 1e8*randn(n); A4 = 1e7*randn(n); so as to have large variation in the norms of its coefficient matrices, the latter being fairly well conditioned (see Table 1).It follows from this table that the set of abscissae of the Newton polygon associated with t × p(x) is K = {0, 1, 3, 4}.Thus t × p(x) has three tropical roots, of multiplicity one, two, and one, respectively.The eigenvalues of P (λ), which we computed with the MATLAB function polyeig, are located in three separate annuli A i , i = 1, 2, 3, given in the first rows of Table 2.These are to be compared to the annuli from the Bini et al. generalized Pellet's theorem (see Theorem 3.1), Melman's generalized Pellet's theorem (see [17,Thm. 3.3] or Theorem 3.1 with q k (x) in place of q k (x)), and that of Theorem 2.7 referred to as Pellet 1, Pellet 2, and Tropical, respectively, in Table 2.The generalized Pellet theorem identifies more annuli with q k (x) in (3.1) than with q k (x) in (3.2), and the bounds provided by the former are tighter as expected from Proposition 3.3.Theorem 2.7 provides only a lower and upper bound for this particular example.It can be seen from Table 2 that the sets H, H, and K, which are defined in (3.3), (3.4), and (2.17), are H = {0, 1, 3, 4}, H = {0, 3, 4}, and K = {0, 4} so that K ⊂ H ⊂ H ⊆ K in accordance with Theorem 3.5.Note that so, for this example, the tropical roots offer an order of magnitude approximation to the eigenvalues of P (λ).Experiment 2. Our next example is a class of matrix polynomials generated via A0 = randn(n); A1 = 1e-3*randn(n); A2 = 1e3*randn(n); A3 = 1e7*randn(n); A4 = 1e-3*randn(n); for a given size n > 3 (but not too large).For this class of matrix polynomials, t × p(x) has only two tropical roots, a small root α 1 of multiplicity three and a large root α 2 of multiplicity one (for n = 5, α 1 = O(10 −3 ) and α 2 = O(10 10 )).Theorems 2.7 and 3.1 detect two annuli, one associated with α 1 containing 3n eigenvalues and one associated with α 2 and containing n eigenvalues.The function polyeig of MATLAB does not find any eigenvalue in the largest annuli.Instead, it tends to return n eigenvalues at infinity.The leading coefficient A 4 is however generically nonsingular, so that there should be no eigenvalue at infinity.
The function polyeig, which solves the polynomial eigenvalue problem via linearization, is not numerically stable [13], [22].The linearization process used in eigensolvers such as polyeig can also affect the sensitivity of the eigenvalues: a well conditioned eigenvalue for P (λ) may be badly conditioned for the linearization [14].As a result, polyeig can return eigenvalues with no digits of accuracy.We note that there is currently no eigensolver for dense matrix polynomials of degree > 2 with guaranteed backward stability.With the aim of addressing this issue, Gaubert and Sharify [9] propose to solve q tropically scaled polynomial eigenvalue problems with the matrix polynomials P (μ) in (2.6) scaled with α i , i = 1, . . ., q.We recall below a version of [9, Algorithm 1].

5
Sort the eigenvalues in modulus from small to large.
7 k = k + nm j .8 end Gaubert and Sharify show experimentally that their algorithm tends to compute eigenpairs with smaller backward errors (see (4.1)) than those computed with the classical approach (i.e., without tropical scaling).Sharify and Tisseur [21] show that amongst the eigenpairs returned by Algorithm 4.1, those with eigenvalues of modulus within order one of α i are computed with small backward errors and their condition numbers are not affected by the linearization process.
We note that Algorithm 4.1 is q times more expensive that polyeig, where q ≤ is the number of distinct tropical roots of t × q(x) but it has better numerical stability properties than polyeig and that it delivers more accurate eigenpairs.It is outside the scope of this paper to develop an efficient eigensolver and also to justify the selection of the computed eigenpairs (see lines 5-6 of Algorithm 4.1).
To illustrate the behavior of Algorithm 4.1, we measure the backward error η p (λ, x) for a computed eigenpair (λ, x) of P (λ) with λ finite and nonzero, with the scaled residual [22] (4.1) We consider the backward error to be small if η(λ, x) ≤ ( n)u.To measure the sensitivity of a simple, finite, and nonzero eigenvalue λ of P where x, y are right and left eigenvectors of P with eigenvalue λ and z, w are right and left eigenvectors of L with eigenvalue λ.Ideally, we would like the linearization L of P to be such that κ L (λ) ≈ κ P (λ).The top plot in Figure 2 shows the backward errors for the computed eigenpairs via polyeig and Algorithm 4.1 for P (λ) generated as in Experiment 2 by setting n = 30 and randn('state',0).The bottom plot displays the ratios between the condition number κ L (λ) of λ as an eigenvalue of L and the condition number κ P (λ) of λ as an eigenvalue of P (λ).In our figure, the x-axis is the eigenvalue index and the eigenvalues  The linearization used by polyeig is the reversal of the first companion linearization of the reversal of P (λ) defined by revP (λ) = λ (P (1/λ).The bottom plot shows that, for this example, when no scaling is applied to P (λ), the linearization process increases the eigenvalue condition numbers by a factor 10 7 but when scaling is used such as in Algorithm 4.1, then κ L (λ) ≈ κ P (λ).
Experiment 3. We consider two quartics from the NLEVP collection of nonlinear eigenvalue problems [5], namely, the butterfly problem and the orr − sommerfeld problem.The coefficient matrices of the butterfly problem are generated as follows: c = kron([1e2 1e-2 1e2 1 1e-3],[1 1]); coeffs = nlevp('butterfly',9,c); A0 = coeffs{1}; A1 = coeffs{2}; A2 = coeffs{3}; A3 = coeffs{4}; A4 = coeffs{5}; Both problems have variations in the norms of their coefficient matrices as shown in Table 3.The moduli of the eigenvalues of these matrix polynomials, the tropical roots of t × p(x) as well as the intervals from the generalized Pellet theorem (Theorem 3.1) which contain the moduli of the eigenvalues of P are all plotted in Figure 3.The backward errors for eigenpairs computed with polyeig and Algorithm 4.1 are plotted in Figure 4, and the ratios between the condition number κ L (λ) of λ as an eigenvalue of the linearization L used by the eigensolvers and the condition number κ P (λ) of λ as an eigenvalue of P (λ) are shown in Figure 5.
For the butterfly problem, K = {0, 2, 3, 4} (see (2.5)) and since A 3 is singular, Theorem 2.7 with Lemma 2.3 and the generalized Pellet's theorems identify two annuli, one containing 2n eigenvalues with magnitude around α 1 = ( A 0 2 / A 2 2 ) 1/2 ≈ 5.1 × 10 −1 , and the second and wider annulus containing the remaining 2n eigenvalues and the two tropical roots α 2 = A 2 2 / A 3 2 ≈ 2.4 × 10 2 and α 3 = A 3 2 / A 4 2 ≈ 4.1 × 10 2 .The top plot in Figure 3 shows that the three tropical roots associated with the butterfly problem are good approximations to the magnitude of the eigenvalues.As a consequence of this and the analysis in [21], the eigenpairs computed by Algorithm 4.1 have small backward errors (see top plot in Figure 4) and the linearization    process used by the eigensolver does not increase the eigenvalue condition numbers (see top of Figure 5).We note that polyeig returns eigenpairs with backward errors as large at 10 −10 and the linearization process increases the eigenvalue condition numbers by a factor 10 10 for the 2n largest eigenvalues.
For the orr − sommerfeld problem, Theorem 2.7 and the generalized Pellet's theorems identify only one annulus.The tropical roots do not offer order of magnitude approximations to all the eigenvalues, in particular for the largest ones (see bottom of Figure 3).Nevertheless, Algorithm 4.1 returns eigenpairs with backward errors all less than 10 −13 ≈ 2( n)u, whereas those returned by quadeig can be as large at 10 −4 (see bottom plot in Figure 4).The linearization process used by the eigensolvers increases the eigenvalue condition numbers by a factor at most 10 5 for Algorithm 4.1 and up to 10 15 for polyeig (see bottom plot in Figure 5).

Concluding remarks.
We have identified sufficient conditions under which the tropical roots of t × p(x) = max 0≤i≤ ( A i x i ) are good order of magnitude approximations to the eigenvalues of P (λ) = i=0 λ i A i .These tropical roots are interesting from the numerical point of view since they are cheap to compute and can be used to define a family of eigenvalue parameter scalings for matrix polynomials that can both improve the backward stability of polynomial eigensolvers based on linearizations and help not to increase the eigenvalue condition numbers of the linearized problem ( see section 4).This is confirmed by the analysis in [21].We anticipate that these tropical roots will help in designing a more numerically stable version of polyeig.
a tropical root.A tropical polynomial of degree has tropical roots counting multiplicities.The multiplicity of a finite root α coincides with the variation of the derivative of the map tp at α, lim →0 d tp dx | x+ − d tp dx | x− .The multiplicity of −∞ as a root of tp is given by lim →0 d tp dx | x+ or, equivalently, by inf{j | a j = −∞}.

•
D(r)) and A(a, b) := {λ ∈ C, a ≤ |λ| ≤ b} denotes a closed annulus centered at 0 with a, b such that 0 < a < b.We write • A(a, b) for the interior of A(a, b).Finally, we denote by (2.8) c 2015 SIAM.Published by SIAM under the terms of the Creative Commons 4.0 license Downloaded 02/09/15 to 130.88.16.37.Redistribution subject to CCBY license

Fig. 2 .
Fig. 2. Problem from Experiment 2. Backward errors of computed eigenpairs and ratios between the eigenvalue condition numbers.

Fig. 3 .
Fig. 3. Butterfly and Orr-Sommerfeld problems.The moduli of the eigenvalues are plotted as " * ".The thick red lines define intervals [t h j , s h j ] from the generalized Pellet theorem containing the |λ j |.The blue dash-dotted lines indicate the value of the tropical roots.
It follows from Corollary 2.8 that if A kj−1 and A kj c 2015 SIAM.Published by SIAM under the terms of the Creative Commons 4.0 license Downloaded 02/09/15 to 130.88.16.37.Redistribution subject to CCBY license If A 0 is nonsingular then q 0 (x) has only one real positive root t 0 and P (λ) has no eigenvalue in the open disk has either no real positive root or two real positive roots s k ≤ t k .In the latter case, P (λ) c 2015 SIAM.Published by SIAM under the terms of the Creative Commons 4.0 license Downloaded 02/09/15 to 130.88.16.37.Redistribution subject to CCBY license has kn eigenvalues in the disk D(s k ) and no eigenvalue in the open annulus • A(s k , t k ).(ii) • D(t 0 i) If A h and A k are nonsingular for some h, k such that 1 ≤ h ≤ k ≤ − 1 and p h (x), p k (x) have exactly two real positive roots s h ≤ t h and s k ≤ t k then P (λ) has n(k − h) eigenvalues in the annulus A(t h , s k ).(ii) If A k is nonsingular for some k such that 1 ≤ k ≤ − 1 and p k has exactly two real positive roots s k ≤ t k then P (λ) has exactly nk eigenvalues inside the annulus A(t 0 , s k ).(iii) If A k is nonsingular for some k such that 1 ≤ k ≤ − 1 and p k has exactly two real positive roots s k ≤ t k then P (λ) has exactly n( − k) eigenvalues inside the annulus A(t k , s ).

Table 2
Pellet and tropical localizations of the spectrum Λ(P ) of P (λ) in Experiment 1.

Table 3
Size n, degree , and norm of the coefficient matrices for the butterfly and orr − sommerfeld problems.increasing order of absolute value.Since polyeig wrongly returns 30 eigenvalues at infinity and the backward error in (4.1) and condition numbers in (4.2) are not defined at infinity, η P (λ, x) and κ L (λ)/κ P (λ) are not plotted for these eigenvalues.The top plot shows that none of the eigenpairs returned by polyeig have a small backward error whereas Algorithm 4.1 returns all eigenpairs with a backward error close to the unit roundoff except for the 3n+1 = 91st eigenvalue |λ 91 | = 7.5×10 8 for which η P (λ 91 , x 91 ) = 1.2 × 10 −13 .The matrix polynomial P (λ) has two tropical roots α 1 = 4.5 × 10 −3 with multiplicity 3 and α 2 = 10 10 with multiplicity 1.The largest tropical root α 2 does not quite provide an order of magnitude approximation to λ 91 and tropical scaling with α 2 yields a slightly too large backward error for the eigenpair (λ 91 , x 91 ).