Stable Perturbations of Nonsymmetric Matrices

James V. Burke and Michael L. Overton

June 1991 Final version for Linear Algebra and its Applications

Abstract

A complex matrix is said to be stable if all its eigenvalues have negative real part. Let J be a Jordan block with zeros on the diagonal and ones on the superdiagonal, and consider analytic matrix perturbations of the form A(ϵ) = J + ϵB + O(ϵ^{2}), where ϵ is real and positive. A necessary condition on B for the stability of A(ϵ) on an interval (0, ϵ_{0}), and a sufficient condition on B for the existence of such a family A(ϵ), is (i) Re tr B ≤ 0; (ii) the sum of the elements on the first subdiagonal of B has nonpositive real part and zero imaginary part; (iii) the sum of the elements on each of the other subdiagonals of B is zero. This result is extended to matrices with any number of nonderogatory eigenvalues on the imaginary axis, and to a stability definition based on the spectral radius. A generalized necessary condition, though not a sufficient condition, applies to arbitrary Jordan structures. The proof of our results uses two important techniques: the Puiseux-Newton diagram and the normal form of Arnold. In the nonderogatory case our main results were obtained by Levantovskii in 1980 using a different proof. Practical implications are discussed.

1 Introduction

A matrix is said to be stable if all its eigenvalues lie in the open left half of the complex plane. A problem of great theoretical interest and practical importance is the following: given a matrix A on the boundary of the set of stable matrices, i.e. with no eigenvalue having positive real part and one or more eigenvalues lying on the imaginary axis, what perturbations to the matrix are associated with stability? To be more specific, suppose that A^{(0)} is an n × n Jordan block J_{0}, with zeros on the diagonal and ones on the main superdiagonal, and consider a real perturbation A(ϵ), with A(0) = J_{0}. It is well known, e.g. [24, Section 2.2], that if

A(ϵ) = [ \begin{matrix}0 & 1 \\ & · & · \\ & & · & · \\ & & & · & 1 \\ ϵ & 0 & · & · & 0 \end{matrix} ]

then the multiple zero eigenvalue splits into n eigenvalues

ϵ^{1/n}ω^{h}, h = 1,…, n,  (1)

where ω is the nth principal root of unity exp(2πi/n). Thus if n > 2, at least one of these eigenvalues has positive real part for all ϵ ≠ 0.

In this paper we give a much stronger result. We show, for example, that if A(ϵ) is a real analytic matrix function of the form

A(ϵ) = J_{0} + ϵB + O(ϵ^{2})

then a necessary condition on B for the stability of A(ϵ) on an interval (0, ϵ_{0}), and a sufficient condition on B for the existence of such a family A(ϵ), is

(i) tr B ≤ 0, i.e. the sum of the diagonal elements of B is less than or equal to zero;

(ii) the sum of the elements on the first subdiagonal of B is less than or equal to zero;

(iii) the sum of the elements on each of the other subdiagonals of B is exactly zero.

Note that (iii) excludes (1) for n > 2. In fact, we show that the jth diagonal in the lower triangle (counting the main diagonal as the first) is associated with perturbations of the form

ϵ^{1/j}ω^{h}, h = 1, ⋯, j,

where ω is the jth principal root of unity. This largely explains (i)–(iii), since eigenvalue splitting separated by angles of less than π in the complex plane cannot be tolerated for stability, while the splitting associated with j = 2, namely in opposite directions in the complex plane, is permitted only along the imaginary axis. The case j = 1 corresponds to an eigenvalue shift, permitted only in the negative direction. Our results are established in quite full generality, allowing A^{(0)} and A(ϵ) to be real or complex and requiring knowledge only of the Jordan form of A^{(0)}. In the case of a complex perturbation of J_{0} (taking ϵ to be real), we obtain, instead of (i), that Re tr B ≤ 0, and, instead of (ii), that the sum of the elements on the first subdiagonal has real part less than or equal to zero and imaginary part equal to zero. This is because an imaginary shift in the eigenvalues is permissible, but even a small rotation of a split along the imaginary axis is not.

In the case that A^{(0)} has one or more derogatory multiple eigenvalues, we obtain a generalized necessary condition but no sufficiency condition. For example, in the nondefective (semisimple) case (e.g. A^{(0)} = 0), conditions of the form (ii) and (iii) vanish (since a perturbation of a nondefective eigenvalue is always Lipschitz), and (i) remains an obvious necessary but insufficient condition for stability. In fact, a stronger necessary and sufficient condition is well known for the nondefective case; if A^{(0)} = 0, it is trivially that the maximum real part of the eigenvalues of B be less than or equal to zero. For a more complete discussion of the semisimple case, see [13,21]. To a large extent the present work arose through attempting to answer the questions raised in Section 6 of [21]. The present paper, together with related work in progress [4], not only answers those questions, but establishes much stronger results than anticipated by [21]. For the much more special symmetric or Hermitian case, see [20].

The proof of our results uses two powerful techniques: the Puiseux-Newton diagram for describing roots of polynomials with analytic coefficients, and the normal form of Arnold for perturbation of multiple eigenvalues. These are both explained in the text.

In the nonderogatory case our main results have already been obtained by Levantovskii [15,16,17], who used a different technique to study C^{1} (continuously differentiable) perturbations. Levantovskii also made use of the Arnold normal form, but he did not use the Puiseux-Newton diagram, which we consider to be a key point for the understanding of these results. (Analyticity is required for the use of the diagram.) See Levantovskii’s papers, as well as [2, p.255-257], for some further results, particularly illustrations of the shape of the stability region in parameter space for various cases.

We note that the Puiseux-Newton diagram has been used in a related context, namely the stability of methods for solving ordinary differential equations; see [11]. However, that work has a very different emphasis, not being concerned with multiple eigenvalues but instead with second-order effects of the perturbation of simple eigenvalues.

All the results have a straightforward extension to the case of reducing the spectral radius of a matrix with multiple eigenvalues. This case is relevant, for example, when the definition of matrix stability requires the magnitude of the eigenvalues to be less than one. This situation typically arises in practice when studying the solution of difference equations rather than differential equations.

The paper is organized as follows. Section 2 uses the Puiseux-Newton diagram to study the roots of polynomials with complex coefficients. Section 3 uses the Arnold normal form to extend these results to perturbation of a complex matrix with multiple eigenvalues, both the nonderogatory and derogatory cases. Section 4 briefly discusses the real case. Section 5 is concerned with reducing the spectral radius of a matrix. Section 6 discusses practical implications.

2 Roots of Polynomials with Complex Coefficients

A complex function β(ϵ) is said to be analytic, equivalently holomorphic, near ϵ = 0 if it may be expanded in a power series in ϵ, convergent in a neighborhood of ϵ = 0. Consider

P (λ, ϵ) = (λ − λ_{0})^{n} + β_{1}(ϵ)(λ − λ_{0})^{n−1} + ⋯ + β_{n}(ϵ) = 0,  (2)

a polynomial equation in λ with analytic coefficients β_{j}(ϵ), satisfying

β_{j}(0) = 0, j = 1,…, n.  (3)

We may therefore write

β_{j}(ϵ) = β_{j}^{(1)}ϵ + β_{j}^{(2)}ϵ^{2} + · · · .  (4)

We shall restrict ϵ to a nontrivial real interval [0, ϵ_{0}].

It is well known, e.g. [3,13], that the roots of (2) are described by series in fractional powers of ϵ. These series are commonly called Puiseux series, since it was Puiseux [22] who established their convergence; however they were derived formally by Newton two centuries earlier. Newton, of course, was concerned only with real coefficients and real roots, but we shall consider the general complex case first for simplicity. We shall obtain the results we need by making use of a diagram devised by Newton for the purpose of calculating coefficients of Puiseux series [19,23]. Since so many calculation techniques go by the name of Newton, we shall follow [3] in calling this diagram the Puiseux-Newton diagram.

Let β̂_{j} = β_{j}^{(ℓ_{j})} be the first nonzero value in the sequence {β_{j}^{(1)}, β_{j}^{(2)},…}. By definition, ℓ_{j} ≥ 1, j = 1,…, n. If β_{j}(ϵ) is identically zero, take ℓ_{j} = ∞; also, since the coefficient of λ^{n} in P (λ, ϵ) is one, take ℓ_{0} = 0, β̂_{0} = 1. Now plot the values ℓ_{j} versus j, and consider the lower boundary of the convex hull of the points plotted. Let s_{j} be the slope of the line on [j, j + 1] forming part of this boundary, j = 0,…, n − 1. Clearly 1/n ≤ s_{0} ≤ s_{1} ≤ ⋯ ≤ s_{n−1}. Figure 1 shows the diagram for the following example (taken from [3]):

n = 3; λ_{0} = 0; β_{1}(ϵ) = ϵ; β_{2}(ϵ) = −ϵ − ϵ^{2}; β_{3}(ϵ) = ϵ^{2} + 2ϵ^{3}.  (5)

We have ℓ_{0} = 0, ℓ_{1} = 1, ℓ_{2} = 1, ℓ_{3} = 2, and so s_{0} = s_{1} = \frac{1}{2}, s_{2} = 1.

[figure: j; l_{j}; X; X; X; X; 1; 2; 3; 2; 1]

Figure 1

Now consider the following “Ansatz” argument. Suppose a root of (2) is to have the form

λ(ϵ) − λ_{0} = αϵ^{p} + ⋯  (6)

where α is nonzero and p is the smallest power of ϵ in the expansion for this root. Substituting (6) into (2), we need

(α^{n}ϵ^{np} + ⋯) + (β̂_{1}ϵ^{ℓ_{1}} + ⋯)(α^{n−1}ϵ^{(n−1)p} + ⋯) + ⋯
+ (β̂_{n−1}ϵ^{ℓ_{n−1}} + ⋯)(αϵ^{p} + ⋯) + (β̂_{n}ϵ^{ℓ_{n}} + ⋯) = 0.

The terms involving the smallest powers of ϵ are among the terms

α^{n}ϵ^{np}, β̂_{1}α^{n−1}ϵ^{ℓ_{1}+(n−1)p},…, β̂_{n−1}αϵ^{ℓ_{n−1}+p}, β̂_{n}ϵ^{ℓ_{n}}.  (7)

For cancellation to take place, at least two terms with the same smallest power of ϵ must appear. Equivalently, p must equal one or more of the slopes s_{0},…, s_{n−1} defined by the Puiseux-Newton diagram. The following discussion will apply to a particular choice of such p. Define f and d by p = s_{f} = … = s_{f+d−1}, so that the line in the diagram with slope p passes from the point (f, ℓ_{f}) to the point (f +d, ℓ_{f+d}). Cancellation of the coefficients of the terms with the smallest powers of ϵ in (7) requires α to be the root of a polynomial equation with degree d, with leading term β̂_{f}α^{d} and constant coefficient β̂_{f+d}, and with an additional intermediate nonzero term for each point (j, ℓ_{j}) lying on the line in the diagram with slope p, where f < j < f + d. Now let p = q/r, where q, r are relatively prime integers. By definition, p is an integral multiple of 1/d, so d is an integral multiple of r, say d = mr. It is then clear from the diagram that of the d − 1 abscissa values j between f and f + d, only every rth value is a candidate for the intersection of the line with a point with integer coordinates. Consequently the polynomial of degree d in α reduces to a polynomial of degree m in α^{r}, which we may denote by Q(γ). The conclusion is that the given value of p is associated with d roots with an expansion of the form (6), with α taking the values

γ_{h}^{1/r}ω^{j}, h = 1,…, m, j = 1,… r  (8)

where the γ_{h} are the m roots of Q(γ) = 0, γ_{h}^{1/r} is the principal rth root of γ_{h} and ω is the principal rth root of unity.

Completing the example given above, we see that the two values for p are s_{0} = s_{1} = \frac{1}{2} and s_{2} = 1. In the case p = \frac{1}{2} we have f = 0, d = 2, r = 2, m = 1, with Q(γ) = γ − 1, so the possible values for α are ±1, giving the Puiseux series

λ(ϵ) − λ_{0} = ±ϵ^{\frac{1}{2}} + · · · .

In the case p = 1 we have f = 2, d = 1, r = 1, m = 1, with Q(γ) = γ − 1, so the only possible value for α is 1, giving the Puiseux series

λ(ϵ) − λ_{0} = ϵ + · · · .

The subsequent terms in the series may also be calculated by repeating the process.

The Puiseux-Newton diagram may be used to establish many results. For example, it is a trivial consequence of the diagram that there is a Puiseux series (6) with p = 1/n if and only if β_{n}^{(1)} is nonzero.

We may now state the main result of this section.

Theorem 1 Consider the polynomial equation (2), with roots given by one or more Puiseux series of the form (6). Suppose that there exists ϵ_{0} > 0 such that all the roots λ(ϵ) of (2) satisfy

Re (λ(ϵ) − λ_{0}) ≤ 0  (9)

for 0 ≤ ϵ ≤ ϵ_{0}. Then

Re β_{1}^{(1)} ≥ 0,  (10)
Re β_{2}^{(1)} ≥ 0, Im β_{2}^{(1)} = 0,  (11)
β_{j}^{(1)} = 0, j = 3,…, n.  (12)

Proof. The coefficient β_{1}(ϵ) is the sum of the differences λ_{0} − λ(ϵ) over the roots λ(ϵ) of (2); thus (10) follows from (9), letting ϵ → 0. The other results follow from the Puiseux-Newton diagram as follows. Consider the Puiseux series corresponding to p = s_{0}, the smallest possible value. In order for the roots on the right-hand side of (6) to all be in the left half-plane, it is necessary, because of (8), that r = 1 or r = 2, with Re γ_{h} ≤ 0, h = 1,…, m in the former case and

Re γ_{h} ≤ 0, Im γ_{h} = 0, h = 1,…, m  (13)

in the latter case. In both cases p ≥ \frac{1}{2}. Since the points (j, ℓ_{j}) lie on or above the line through the origin with slope p = s_{0}, (12) holds. In the case r = 1, p ≥ 1 so β_{2}^{(1)} is also zero. In the case r = 2, we have

Q(γ) = γ^{m} + β_{2}^{(1)}γ^{m−1} + ⋯ + β_{d}^{(m)}.  (14)

Since −β_{2}^{(1)} is the sum of the roots of Q(γ), (11) follows from (13). □

The result of Theorem 1 was first obtained by Levantovskii, using a completely different technique developed for C^{1} polynomial coefficients. We were unable to verify some of the steps in his proof, particularly the justification for the division by the quantity in [15, p.20, line D]. Although we do not doubt that the result is correct for C^{1} perturbations, we think that our approach using the Puiseux-Newton diagram, which requires analyticity, gives far more insight.

The next theorem shows that the conditions (10)–(12), as well as being necessary for (9), are sufficient for the existence of a polynomial whose roots have real part strictly less than Re λ_{0}.

Theorem 2 Suppose that (10)–(12) hold. Then there exists a polynomial of the form (2) such that

Re (λ(ϵ) − λ_{0}) < 0

on an open interval (0, ϵ_{0}).

Proof. The proof is trivial if n = 1. Otherwise, let β_{k}(ϵ), k = 1,…, n, be defined by the coefficients of

(λ−λ_{0} + \frac{β_{1}^{(1)}}{n} ϵ+ϵ^{2})^{n−2}(λ−λ_{0} +i(β_{2}^{(1)}ϵ)^{\frac{1}{2}} + \frac{β_{1}^{(1)}}{n} ϵ+ϵ^{2})(λ−λ_{0} −i(β_{2}^{(1)}ϵ)^{\frac{1}{2}} + \frac{β_{1}^{(1)}}{n} ϵ+ϵ^{2}).

(Here i = \sqrt{−1}.) □

Note that two roots split from the rest, but that all the roots have the same real part. (It is interesting to compare this with the more complicated example in [15, p.20, line F], as well as with [21, p. 489].)

It is also possible to use the Puiseux-Newton diagram to make statements about higher order coefficients of β_{j}(ϵ). For example, it seems reasonable to suppose that the following theorem was known to Newton, since its proof using the diagram is trivial.

Theorem 3 Let g be the minimum, over all the Puiseux series describing the roots, of the associated values p. Then

β_{j}^{(k)} = 0, for all j > ⌊k/g⌋.  (15)

Here ⌊x⌋ is the largest integer less than or equal to x. This theorem may be viewed as a form of converse to the obvious fact that, if β_{j}(ϵ) is identically zero, j = n − k,…, n, then for any of the Puiseux series describing the roots, p ≥ 1/k, since dividing P (λ, ϵ) by the factor λ^{n−k} leaves a polynomial of degree k.

Finally we note that Theorem 1 can be generalized in a straightforward way to products of polynomials of the form (2). Since our main interest is in matrices, we omit the details (but see [4]).

3 Eigenvalues of Complex Matrices

Let A^{(0)} ∈ C^{n×n} be a matrix with eigenvalues λ_{1},…, λ_{σ}, with (algebraic) multiplicities t_{1},…, t_{σ} respectively, where Re λ_{1} = … = Re λ_{ρ} = 0 and Re λ_{ρ+1} < 0,…, Re λ_{σ} < 0. Let the Jordan form of A^{(0)} be given by

A^{(0)} = P JP^{−1}

where

J = [ \begin{matrix}J_{1} \\ & . & . & . \\ & & & & J_{σ} \end{matrix} ] ,

J_{k} = [ \begin{matrix}J_{k1} \\ & . & . & . \\ & & & & J_{km_{k}} \end{matrix} ] ,

and the Jordan block

J_{kl} = [ \begin{matrix}λ_{k} & 1 \\ & · & · \\ & & · & · \\ & & & · & 1 \\ & & & & λ_{k} \end{matrix} ]

has dimension n_{kl}. We have

n_{k1} + ⋯ + n_{km_{k}} = t_{k}, k = 1,…, σ.

If m_{k} = 1, λ_{k} is said to be a nonderogatory eigenvalue, while if m_{k} = t_{k}, i.e. n_{k1} = ⋯ = n_{km_{k}} = 1, λ_{k} is said to be nondefective (semisimple).

We will now consider an analytic perturbation of the matrix A^{(0)}, i.e. a matrix A(ϵ) ∈ C^{n×n} each of whose elements is an analytic function of ϵ near ϵ = 0, with A(0) = A^{(0)}. It will be convenient to define

J(ϵ) = P^{−1}A(ϵ)P

so that J(0) = J, the Jordan form of A^{(0)}.

We will use a powerful result of V.I. Arnold [1, Section 1 together with Theorem 4.4] which states that any such analytic perturbation J(ϵ) of J has the following local representation:

J(ϵ) = Y (ϵ)C(ϵ)Y (ϵ)^{−1}

where Y (ϵ) is analytic, with Y (0) = I,

C(ϵ) = Diag(C_{k}(ϵ)),
C_{k}(ϵ) = J_{k} + E_{k}(ϵ),

and E_{k}(ϵ), which satisfies E_{k}(0) = 0, has a structure best described by an example as follows:

[ \begin{matrix}X & X & X & X & X & X \\ X \\ X & & & X & X & X \\ X & & & X & & X \end{matrix} ] .

Here we suppose that J_{k} has three blocks J_{k1}, J_{k2}, J_{k3}, with dimensions n_{k1} = 3, n_{k2} = 2, n_{k3} = 1 respectively, so that t_{k} = 6. Each X indicates a distinct complex analytic function of ϵ. Denote the diagonal blocks of C_{k}(ϵ) by C_{kl}(ϵ); each of these has the form commonly known as a companion matrix. Let us write

C_{kl}(ϵ) = [ \begin{matrix}λ_{k} & 1 \\ & · & · \\ & & · & · \\ & & & · & · \\ & & & & λ_{k} & 1 \\ −θ_{kl}^{(m_{k})}(ϵ) & · & · & · & −θ_{kl}^{(2)}(ϵ) & λ_{k} − θ_{kl}^{(1)}(ϵ) \end{matrix} ] ,  (16)

so that, for the example given above,

C_{k}(ϵ) = [ \begin{matrix}λ_{k} & 1 \\ & λ_{k} & 1 \\ −θ_{k1}^{(3)}(ϵ) & −θ_{k1}^{(2)}(ϵ) & λ_{k} − θ_{k1}^{(1)}(ϵ) & X & X & X \\ X & & & λ_{k} & 1 \\ X & & & −θ_{k2}^{(2)}(ϵ) & λ_{k} − θ_{k2}^{(1)}(ϵ) & X \\ X & & & X & & λ_{k} − θ_{k3}^{(1)}(ϵ) \end{matrix} ] .

Again each X indicates a distinct unnamed analytic function of ϵ. Note C_{k}(0) = J_{k}.

The matrix C(ϵ) is called a versal deformation of J, and the proof that this versal deformation exists relies on some elementary concepts from differential geometry, in particular the orbit of a matrix and the notion of transversality with respect to the orbit. The orbit of J is simply the set of matrices which are similar to J, i.e. have exactly the same Jordan form. The orbit is a differential manifold in the matrix space C^{n×n} and so there exists a linear manifold T which is tangent to the orbit at the point J. The main idea is as follows: since J is similar to all the matrices in its orbit, any perturbation J(ϵ) is similar to a representation C(ϵ) which lies entirely in a linear manifold of matrices transversal to the orbit at J, e.g. the orthogonal complement of T. For example, consider the “least generic” case, where J is the identity matrix (which has exactly one nondefective eigenvalue); then the orbit consists of the single point I and the orthogonal complement of T is C^{n×n}. Thus C(ϵ) is a matrix with n^{2} distinct analytic component functions. On the other hand, the most generic case is where J is a diagonal matrix with distinct elements λ_{1},…, λ_{n}. In this case it turns out that the orthogonal complement to T is the set of all diagonal matrices, so that it suffices for C(ϵ) to contain only n distinct analytic component functions, the offdiagonal elements being zero. In general it turns out that the orthogonal complement to T is the set of matrices which commute with the transpose of J, which leads to the general minimum-parameter form described above. (This particular normal form is not in the orthogonal complement, but in another transversal linear manifold.) For further details, see either Arnold’s original paper, or the excellent summary given by Fairgrieve [6].

Before stating the main results we need a definition.

Definition. Define the jth generalized trace of a square matrix A, denoted by

tr^{(j)}A,

as the sum of the elements on that diagonal of A which is j − 1 positions below the main diagonal. Thus one obtains the ordinary trace in the case j = 1 and the bottom left element of the matrix in the case that j is the dimension of the matrix. If j exceeds the dimension of A, take tr^{(j)}A = 0.

Theorem 4 Let A(ϵ) be an analytic matrix function with A(0) = A^{(0)}, A′(0) = A^{(1)}. We shall restrict ϵ to a nontrivial real interval [0, ϵ_{0}]. Partition P^{−1}A^{(1)}P conformally with the partition of J given above and denote its diagonal block corresponding to J_{k} by B_{k}, k = 1,…, σ, with each B_{k} having diagonal blocks B_{kl} corresponding to J_{kl}, l = 1,…, m_{k}. Define

τ_{k}^{(j)} = ∑_{l=1}^{m_{k}} tr^{(j)}B_{kl}, j = 1,…, max_{l}{n_{kl}}, k = 1,…, σ.

(This sum includes zero terms when j exceeds n_{kl}.) Then, for k = 1,…, σ, the eigenvalues of A(ϵ) corresponding to λ_{k} are the roots of

(λ_{k} − λ_{0})^{t_{k}} + β_{k1}(ϵ)(λ_{k} − λ_{0})^{t_{k}−1} + ⋯ + β_{kt_{k}}(ϵ) = 0  (17)

where the β_{kj}(ϵ) are analytic, with

β′_{kj}(0) = −τ_{k}^{(j)}, j = 1,…, max_{l}{n_{kl}}.

Proof. Define J(ϵ) = P^{−1}A(ϵ)P, and let C(ϵ) be the versal deformation of J(ϵ) described above. We have

J′(0) = C′(0) + Y′(0)J + J(Y^{−1})′(0)
= C′(0) + [Y′(0), J],  (18)

where the latter term is the Lie bracket

[Y′(0), J] = Y′(0)J − JY′(0).

Following Arnold, we note that

⟨[Y′(0), J], K⟩ = 0  (19)

for all K such that K^{∗} commutes with J, where

⟨A, B⟩ = tr A^{∗}B,

the Frobenius matrix inner product, and A^{∗} denotes the complex conjugate transpose of A. Such matrices K^{∗} are block diagonal, K^{∗} = Diag(K_{k}^{∗}), where K_{k}^{∗} commutes with J_{k}. From [1,8], the matrices commuting with J_{k} have a block upper triangular Toeplitz form, specifically, for the example given above,

K_{k}^{∗} = [ \begin{matrix}a & b & c & d & e & f \\ & a & b & & d \\ & & a \\ & g & h & i & j & k \\ & & g & & i \\ & & l & & m & n \end{matrix} ] .

Consequently, K_{k} has a block lower triangular Toeplitz structure. Note that if m_{k} = 1 (the nonderogatory case), K_{k} is a lower triangular Toeplitz matrix, while if m_{k} = t_{k} (the nondefective case), there is no restriction on the structure of K_{k}. By (19), then, we have

tr^{(j)}G_{kl} = 0, j = 1,…, max_{l}{n_{kl}}, k = 1,…, σ,

where the G_{kl} are the diagonal blocks of [Y′(0), J], corresponding in dimension to B_{kl} and J_{kl}. (Indeed, more can be said about the off-diagonal blocks of [Y′(0), J] but we shall not need this.) Now, by definition, the diagonal blocks of J′(0) are B_{kl}, so, by (18),

τ_{k}^{(j)} = ∑_{l=1}^{m_{k}} tr^{(j)}C′_{kl}(0) = − ∑_{l=1}^{m_{k}} θ_{kl}^{(j)′}(0);  (20)

here θ_{kl}^{(j)} is taken to be zero if j exceeds n_{kl}.

Since A(ϵ) is similar to J(ϵ) and therefore to C(ϵ), we need only examine the eigenvalues of C(ϵ). We need the following:

Lemma 1

det (λI − C_{k}(ϵ)) = (λ − λ_{k})^{t_{k}} + β_{k1}(ϵ)(λ − λ_{k})^{t_{k}−1} + … + β_{kt_{k}}(ϵ)  (21)

where

β′_{kj}(0) = ∑_{l=1}^{m_{k}} θ_{kl}^{(j)′}(0).

In the case that λ_{k} is nonderogatory, the proof of Lemma 1 is trivial. In general, it may be proved by a fairly straightforward but complicated induction argument which gives little insight; we shall not give it here. For a more direct but perhaps less intuitive proof, see [4]. The proof of Theorem 4 is completed by combining (20) and (21). □

The next theorem continues to use the notation of Theorem 4.

Theorem 5 Suppose now that all the eigenvalues of A(ϵ) have nonpositive real part on a nontrivial interval [0, ϵ_{0}]. Then, for k = 1,…, ρ,

Re τ_{k}^{(1)} ≤ 0,  (22)

Re τ_{k}^{(2)} ≤ 0, Im τ_{k}^{(2)} = 0,  (23)

and

τ_{k}^{(j)} = 0, j = 3,…, max_{l}{n_{kl}}.  (24)

Proof. First note that there is no restriction on the eigenvalues of C_{k}(ϵ) for k > ρ as long as ϵ_{0} is small enough, since Re λ_{k} < 0 for k > ρ. The proof is completed by applying Theorem 1 to equation (17), k = 1,…, ρ. □

If λ_{k} is nonderogatory, i.e. m_{k} = 1 and J_{k} = J_{k1}, B_{k} = B_{k1}, we have τ_{k}^{(j)} = tr^{(j)}B_{k}, j = 1,…, t_{k}. The next theorem shows that the necessary condition of Theorem 5 is also a sufficient condition for the existence of a stable family A(ϵ) when the eigenvalues with real part equal to zero are all nonderogatory.

Theorem 6 Assume that λ_{k} is nonderogatory, i.e. m_{k} = 1, for k = 1,…, ρ. Let B be a given matrix and denote its diagonal blocks, partitioned conformally with J, by B_{k}, k = 1,…, σ. Suppose that B satisfies

Re tr ^{(1)}B_{k} ≤ 0,
Re tr ^{(2)}B_{k} ≤ 0, Im tr ^{(2)}B_{k} = 0,
tr ^{(j)}B_{k} = 0, j = 3,…, t_{k},

for k = 1,…, ρ. Then there exists an analytic family A(ϵ), with A(0) = A^{(0)}, A′(0) = P BP^{−1}, such that the maximum real part of the eigenvalues of A(ϵ) is strictly negative on an open interval (0, ϵ_{0}).

Proof. There exists an Arnold normal form for J + ϵB, namely

J + ϵB = Y (ϵ)C(ϵ)Y (ϵ)^{−1},

where Y (ϵ) is analytic and satisfies Y (0) = I, and where C(ϵ) = Diag(C_{k}(ϵ)). We need consider only k ≤ ρ as long as ϵ is small enough. By the nonderogatory assumption, C_{k}(ϵ) consists of a single companion matrix block (see (16)). We have

B = C′(0) + [Y′(0), J].

Observe that the relevant diagonal blocks of the second term have generalized traces equal to zero, so the first-order terms in the last row of C_{k}(ϵ) are tr^{(j)}B_{k}, j = t_{k},…, 1. Now define C̃(ϵ) to be another block diagonal matrix with the same structure as C(ϵ). Specifically, choose C̃(ϵ) so that C̃(0) = C(0) = J and C̃′(0) = C′(0), with the kth block of C̃(ϵ) set equal to the companion matrix for the corresponding stable polynomial (of degree t_{k}) defined by Theorem 2. Let

J(ϵ) = J + ϵB + Y (ϵ)(C̃(ϵ) − C(ϵ))Y (ϵ)^{−1}
= Y (ϵ)C̃(ϵ)Y (ϵ)^{−1}

so that J′(0) = B and, by Theorem 2, J(ϵ) is stable on an interval (0, ϵ_{0}). The proof is completed by defining

A(ϵ) = P J(ϵ)P^{−1}. □

We note that, in the nonderogatory case, Theorem 4 is implicit in the work of Levantovskii. However, he apparently did not obtain our necessary condition in the derogatory case. Regarding analyticity, Arnold [1, p.37] mentions that his normal form may be extended to C^{1} perturbations, although the proof is tedious. As in the polynomial case, this extension does not seem to be worth the effort involved.

Having discussed the nonderogatory case, consider now the other extreme case, namely when all the eigenvalues with real part equal to zero are nondefective (semisimple), i.e. m_{k} = t_{k}, k = 1,…, ρ. The necessary condition of Theorem 5 reduces to

Re tr B_{k} ≤ 0, k = 1,…, ρ

since all the B_{kl} are of dimension one. This condition is an obvious necessary condition for stability, and equally clearly not sufficient for the existence of a stable family A(ϵ) with A′(0) = P BP^{−1}. A necessary and sufficient condition is, however, well known in this case, as the following theorem states.

Theorem 7 Assume that λ_{k} is nondefective, i.e. m_{k} = t_{k}, for k = 1,…, ρ. Let A(ϵ) be any matrix family with

A(0) = A^{(0)}, A′(0) = P BP^{−1},  (25)

and where, as before, B_{k} is the diagonal block of B corresponding to λ_{k}. Then a necessary condition for all the eigenvalues of A(ϵ) to have nonpositive real part on a nontrivial real interval [0, ϵ_{0}] [respectively sufficient condition for all the eigenvalues of A(ϵ) to have negative real part on an open interval (0, ϵ_{0})] is that, for k = 1,…, ρ, the maximum real part of the eigenvalues of B_{k} is less than or equal to zero [respectively less than zero].

The proof uses [13,21]. Note that in this case the relevant first ρ columns of P are eigenvectors. Note also that, unlike in Theorem 6, the sufficient condition applies to any family A(ϵ) satisfying (25); this is because the eigenvalues are Lipschitz in ϵ. Analyticity of A(ϵ) is not required; C^{1} parameter dependence is sufficient.

When λ_{k} is simple, i.e. m_{k} = t_{k} = 1, k = 1,…, ρ, the necessary conditions of Theorems 5 and 7 and the sufficient condition of Theorem 6 all reduce to the following: the first ρ diagonal elements of B have nonpositive real part.

The most difficult case is therefore the derogatory, defective case, for which we have given only a necessary condition. Perhaps a sufficient condition could be derived by consideration of the off-diagonal blocks of [Y′(0), J].

4 The Real Case

Every eigenvalue of a real matrix is either real or one of a complex conjugate pair. The results given in Section 3 apply, in particular, in the real case. However, they may be simplified. When A^{(0)} is real, its Jordan form may be taken to have the form

A^{(0)} = P JQ^{T}

where P^{T} Q = I and

J = Diag(J_{0}, J_{1}, \overline{J}_{1},…, J_{s}, \overline{J}_{s},…);
P = [P_{0}, P_{1}, \overline{P}_{1},…, P_{s}, \overline{P}_{s},…];
Q = [Q_{0}, Q_{1}, \overline{Q}_{1},…, Q_{s}, \overline{Q}_{s},…].

Here J_{0} corresponds to the eigenvalue λ_{0} = 0, the only possible real eigenvalue on the imaginary axis; this block and the corresponding blocks P_{0}, Q_{0} may be absent. The blocks P_{0} and Q_{0} can be taken to be real. Blocks J_{1},…, J_{s} correspond to eigenvalues with positive imaginary part and zero real part. The symbol \overline{z}, of course, denotes the complex conjugate of z. The number of eigenvalues on the imaginary axis is ρ = 2s + 1 if the zero eigenvalue is present and ρ = 2s otherwise. Blocks not explicitly listed are of no interest since their eigenvalues have negative real part.

Now assume also that A(ϵ) is real. We have, if the zero eigenvalue is present, that B_{0} = Q_{0}^{T} A^{(1)}P_{0} is real, and hence that τ_{0}^{(j)} is real for all j. Consequently, conditions (22)–(24) reduce, in the case k = 0, to τ_{0}^{(1)} ≤ 0, τ_{0}^{(2)} ≤ 0, and τ_{0}^{(j)} = 0 for j > 2.

Now consider a complex conjugate pair of eigenvalues. In this case the τ_{k}^{(j)} are generally complex, but they occur in complex conjugate pairs. Since each of conditions (22)–(24) applies equally to a complex number or its conjugate, each condition appears exactly twice. Thus it is necessary to apply them only to τ_{k}^{(j)} corresponding to eigenvalues with positive imaginary part.

Most likely it is possible to rewrite these conditions using only real quantities by using the real Jordan form [7,14]. However, our attempts to do this have led only to complicated results. See [7] for the extension of the Arnold normal form to the real field.

Regarding complex conjugate roots of polynomials with real coefficients: instead of (2) we may consider

P (λ, ϵ) = [(λ − ia)^{n} + β_{1}(ϵ)(λ − ia)^{n−1} + ⋯ + β_{n}(ϵ)] ·
[(λ + ia)^{n} + β̄_{1}(ϵ)(λ + ia)^{n−1} + ⋯ + β̄_{n}(ϵ)]  (26)

where a is real and nonzero. Suppose now that P (λ, ϵ) has no root with positive real part on a nontrivial interval [0, ϵ_{0}]. By the generalization of Theorem 1 mentioned at the end of Section 2, it follows that (10)–(12) hold. Let

R(λ) = λ^{2} + a^{2}.

Multiplying out the factors in (26) we see that

P (λ, ϵ) = R(λ)^{n−2}(R(λ)^{2} + ν_{1}ϵλR(λ) + ν_{2}ϵR(λ) + ν_{3}ϵ) + O(ϵ^{2})

where ν_{1} ≥ 0, ν_{3} ≤ 0, and there is no sign restriction on ν_{2}. The same result (with a different derivation) is given in [15, p.20, line E].

5 Reducing the Spectral Radius

The spectral radius of a matrix is the radius of the smallest disk in the complex plane containing all its eigenvalues. For the purposes of this section, change the definition of matrix stability to require the spectral radius to be less than one. We first modify the polynomial theorems accordingly.

Theorem 8 Consider the polynomial equation (2), with roots given by one or more Puiseux series of the form (6). Suppose that there exists ϵ_{0} > 0 such that all the roots λ(ϵ) of (2) satisfy

|λ(ϵ)| ≤ |λ_{0}|  (27)

for 0 ≤ ϵ ≤ ϵ_{0}. Then

Re \overline{λ}_{0}β_{1}^{(1)} ≥ 0,  (28)
Re \overline{λ}_{0}^{2}β_{2}^{(1)} ≥ 0, Im \overline{λ}_{0}^{2}β_{2}^{(1)} = 0,  (29)
β_{j}^{(1)} = 0, j = 3,…, n.  (30)

Proof. First note that when λ_{0} is real and positive the result is the same as before; this is because the circle of radius |λ_{0}| centered at the origin lies to the left of the vertical axis Re z = Re λ_{0}. When λ_{0} is real and negative we obtain the negative of the original conditions, since the circle lies to the right of the corresponding axis. In general, (27) is equivalent to

Re \overline{λ}_{0}(λ(ϵ) − λ_{0}) ≤ −\frac{1}{2} |λ(ϵ) − λ_{0}| ^{2}

so that a necessary condition for (27) to hold is that

Re \overline{λ}_{0}(λ(ϵ) − λ_{0}) ≤ 0.  (31)

The proof is now an easy modification of the proof of Theorem 1 as follows. The coefficient β_{1}(ϵ) is the sum of the differences λ_{0} − λ(ϵ) over the roots λ(ϵ) of (2); thus the first condition follows from (27), letting ϵ → 0. The other results follow from the Puiseux-Newton diagram. In order for the roots on the right-hand side of (6) to lie in the half-plane defined by (31), it is necessary that r = 1 or r = 2, with Re \overline{λ}_{0}γ_{h} ≤ 0 in the former case and

Re \overline{λ}_{0}^{2}γ_{h} ≤ 0, Im \overline{λ}_{0}^{2}γ_{h} = 0, h = 1,…, m

in the latter case. The rest of the proof is exactly as before. □

A sufficient condition is somewhat harder to obtain in the spectral radius case. The difficulty is, by way of example, that an eigenvalue splitting of the form 1 ± iϵ^{\frac{1}{2}} , while having the same real part as the unperturbed eigenvalue 1, has spectral radius \sqrt{1 + ϵ}. Consequently, we need to make an assumption that the left-hand side of (28) is sufficiently bigger than the first left-hand side of (29). For example, we have:

Theorem 9 Suppose that λ_{0} ≠ 0 and that (28)–(30) hold. Suppose also that

|β_{2}^{(1)}| < 2Re \overline{λ}_{0} \frac{β_{1}^{(1)}}{n} .  (32)

Then there exists a polynomial of the form (2) such that

|λ(ϵ)| < |λ_{0}|

on an open interval (0, ϵ_{0}).

Proof. Let β_{k}(ϵ), k = 1,…, n, be defined by the coefficients of

(λ − λ_{0} + \frac{β_{1}^{(1)}}{n} ϵ)^{n−2}(λ − λ_{0} + i(β_{2}^{(1)}ϵ)^{\frac{1}{2}} + \frac{β_{1}^{(1)}}{n} ϵ)(λ − λ_{0} − i(β_{2}^{(1)}ϵ)^{\frac{1}{2}} + \frac{β_{1}^{(1)}}{n} ϵ). □

Note that the larger |λ_{0}| is, the more likely (32) is satisfied for given β_{1}^{(1)}, β_{2}^{(1)}. This reflects the fact that the the smaller a circle is, the greater the distance between it and a tangent line segment of given length. In the limit case λ_{0} = 0, reduction of the spectral radius is, of course, impossible. The factor n in (32) reflects the particular choice of polynomial used in Theorem 2 and can be improved (see [4]).

Now take A^{(0)} to have the same Jordan form as in Section 3, but assuming |λ_{1}| = … = |λ_{ρ}| = 1 and |λ_{ρ+1}| < 1,…, |λ_{σ}| < 1. Instead of Theorem 5, we obtain:

Theorem 10 Suppose that the spectral radius of A(ϵ) is less than or equal to one on a nontrivial interval [0, ϵ_{0}]. Then, for k = 1,…, ρ,

Re \overline{λ}_{k}τ_{k}^{(1)} ≤ 0,
Re \overline{λ}_{k}^{2}τ_{k}^{(2)} ≤ 0, Im \overline{λ}_{k}^{2}τ_{k}^{(2)} = 0,

and

τ_{k}^{(j)} = 0, j = 3,…, max_{l}{n_{kl}}.

The proof is identical to the proof of Theorem 5, using Theorem 8 instead of Theorem 1. In the nonderogatory case, we obtain the following sufficient condition.

Theorem 11 Assume that λ_{k} is nonderogatory, i.e. m_{k} = 1, for k = 1,…, ρ. Let B be a given matrix and denote its diagonal blocks, partitioned conformally with J, by B_{k}, k = 1,…, σ. Suppose that B satisfies

Re \overline{λ}_{k}tr ^{(1)}B_{k} ≤ 0,
Re \overline{λ}_{k}^{2}tr ^{(2)}B_{k} ≤ 0, Im \overline{λ}_{k}^{2}tr ^{(2)}B_{k} = 0,
tr ^{(j)}B_{k} = 0, j = 3,…, t_{k},

for k = 1,…, ρ. Suppose also that

|tr ^{(2)}B_{k}| < −\frac{2}{n}Re \overline{λ}_{k}tr ^{(1)}B_{k}.  (33)

Then there exists an analytic family A(ϵ), with A(0) = A^{(0)}, A′(0) = P BP^{−1}, such that the spectral radius of A(ϵ) is strictly less than one on an open interval (0, ϵ_{0}).

The proof combines those of Theorem 6 and 9. As before, (33) can be weakened. Finally, in the nondefective case we have the following modification to Theorem 7.

Theorem 12 A necessary condition for the spectral radius of A(ϵ) to be less than or equal to one on a nontrivial real interval [0, ϵ_{0}] [respectively sufficient condition for the spectral radius of A(ϵ) to be less than one on an open interval (0, ϵ_{0})] is that, for k = 1,…, ρ, the maximum real part of the eigenvalues of \overline{λ}_{k}B_{k} is less than or equal to zero [respectively less than zero].

The proof uses [21].

6 Practical Implications

How useful are our results in practice? The necessary condition is practical in the sense that, if the Jordan form of A^{(0)} is known, a candidate first-order perturbation A^{(1)} can be rejected immediately if it does not satisfy the given necessary condition. However, the sufficient condition is of little practical use, since nothing is said about how to compute the actual matrix family A(ϵ). Indeed, it is well known that the nonsymmetric eigenvalue problem is very ill-conditioned in the case of multiple eigenvalues. The reason for the ill-conditioning is immediately apparent from the Puiseux series (6); the eigenvalues are not Lipschitz when the Puiseux exponent is less than one. Computing the Jordan form of A^{(0)} is itself a very ill-conditioned numerical problem and, in fact, it is the instability of this process which originally led to Arnold’s normal form for perturbed matrices.

Instead of the Jordan form, one may consider reduction of A^{(0)} to some other canonical form. The most stable of these is the triangular or Schur form, since only unitary similarity transformations are required, but the Schur form gives little information beyond the eigenvalues themselves. Nonunitary transformations are required to further reduce a matrix to block diagonal form, with one triangular block associated with each eigenvalue [10, Section 7.1.3]. In fact, our necessary condition (but apparently not the sufficient condition) can be extended to apply to the block diagonal canonical form [4].

The idea of the Schur or block diagonal forms, as opposed to the Jordan form, is to avoid the reduction to too small a number of parameters. On the other hand the Arnold theory uses minimum-parameter versal deformations, hence the use of the Jordan form and the block triangular Toeplitz matrices which commute with the Jordan form. Several proposals have been made to stabilize the computation of the Jordan form of a matrix A^{(0)}; see [5,9,12]. One approach is to compute the “least generic” of the Jordan forms of those matrices which are within some tolerance of A^{(0)}. Perhaps the most hopeful approach along these lines is that of Fairgrieve [6], who has introduced a path following method for this purpose. Fairgrieve’s method is based on Arnold’s normal form and in particular uses Arnold’s classification of the degeneracy hierarchy [1], in which, as is well known, a nonderogatory eigenvalue is the most generic and a nondefective one the least generic, but for which the ordering of the “middle” cases is by no means obvious.

We conclude with the description of two numerical experiments, conducted in Matlab [18], which give substantial insight. In the first experiment, we set J to the Jordan block of order n = 10 with zeros on the diagonal and ones on the superdiagonal. We then computed the maximum real part of the eigenvalues of J + ϵB^{(k)}, say r_{k}, where ϵ = 10^{−10} and B^{(k)} is a random real matrix of order n with the property that

tr^{(j)}B^{(k)} = 0, j = n − k + 1,…, n.  (34)

Thus, for k = 1, the bottom left element of B^{(k)} is zero; in addition, for k = 2, the two elements on the (n−1)th subdiagonal sum to zero, etc. This matrix was computed using the Matlab pseudo-random number generator, first producing a matrix with elements independently uniformly distributed in the interval [−1, 1], then adding the necessary k constant diagonals to enforce (34). The results were averaged over 100 random matrices and the resulting averaged values of r_{k} plotted, as a function of k, as “∗”’s in Figure 2. The “+” symbols in Figure 2 show the quantities ϵ^{1/(n−k)}. Inspecting Figure 2, one sees a steady reduction in the averaged value of r_{k} as k increases from 0 to n/2, with close agreement between the two plotted values. At k = n/2, however, a “floor” is reached. Variation of ϵ and the machine precision showed that this floor is not a consequence of round-off error, as might be suspected. Instead, it is easily explained by the Puiseux-Newton diagram. By imposing (34), we are setting the corresponding polynomial coefficients β_{j}^{(1)} to zero. Thus for each k between 1 and n/2, slopes with value 1/n, 1/(n − 1),…, 1/(n − k + 1) are not possible. Consequently, roots of the form (6) with corresponding powers are not possible, and the perturbed eigenvalues are instead dominated by terms of order ϵ^{1/(n−k)}. Hence the close association of the “∗” and “+” plots. However, when k reaches n/2, the second-order coefficient β_{n}^{(2)} comes into play. Since the second-order terms arising from the Arnold normal form are not zero, the slope with value 2/n cannot be eliminated from the Puiseux-Newton diagram. Hence the “floor”. We have succeeded only in reducing the order of the eigenvalue perturbations from ϵ^{1/n} to ϵ^{2/n}.

In the second experiment, we set J to a full Jordan block of order n and defined B to be a real random matrix (in the same sense as above) with the property that

tr^{(j)}B = 0, j = 1; j = 3,…, n.

This time we plotted the maximum real part of the eigenvalues of J + ϵB against the value of tr ^{(2)}B. Again, the results were averaged over 100 random cases. Plots are shown in Figure 3 for n = 4 and n = 5. The reason for the two very different graphs is similar to the explanation for the result of the first experiment. In the case n = 5, 19 the nonzero second-order coefficient β_{5}^{(2)} yields a slope with value 2/5 in the Puiseux-Newton diagram; this prevents a reduction in the maximum real part of the perturbed eigenvalues. However, in the case n = 4, the slopes in the diagram must be at least 1/2, so that stability is possible when tr ^{(2)}B ≤ 0. Theorem 5 suggests that the larger the quantity −tr ^{(2)}B is, the “more stable” the perturbed matrix should be, and indeed this is supported by the graph for n = 4. Stability is not actually achieved, however, because of the effect of β_{4}^{(2)} which is, in general, nonzero.

Acknowledgements. We thank V.I. Arnold for bringing the work of Levantovskii to our attention. We also thank R.S. Womersley for some helpful discussions during the early stages of this work. This research was supported in part by National Science Foundation Grants No. DMS-88-03206 and CCR-88-02408.

References

[1] V. I. Arnold, On matrices depending on parameters, Russian Math. Surveys 26:2(1971), 29-43.

[2] V. I. Arnold, Geometrical Methods in the Theory of Ordinary Differential Equations, Springer Verlag, New York and Berlin, (1983).

[3] H. Baumgärtel, Analytic Perturbation Theory for Matrices and Operators, Birkhäuser Verlag, Boston (1985).

[4] J. V. Burke and M. L. Overton, Differential properties of eigenvalues, in preparation.

[5] J.W. Demmel, Computing stable eigendecompositions of matrices, Lin. Alg. Appl. 79(1986), 163-193.

[6] T. F. Fairgrieve, The Application of Singularity Theory to the Computation of the Jordan Canonical Form, M. Sc. Thesis, Department of Computer Science, University of Toronto, Toronto (1986).

[7] D. M. Galin, On real matrices depending on parameters, Uspekhi Math. Nauka 27(1974), 241-242. (In Russian).

[8] F. R. Gantmacher, The Theory of Matrices, volumes 1 and 2, Chelsea, New York (1959).

[9] G. H. Golub and J.H. Wilkinson, Ill-conditioned eigensystems and the computation of the Jordan canonical form, SIAM Review 18 (1976), 578-619.

[10] G.H. Golub and C. Van Loan, Matrix Computations, Johns Hopkins University Press, 2nd edition (1989).

[11] R. Jeltsch, Multistep methods using higher derivatives and damping at infinity, Math. Comp. 31 (1977), 124-138.

[12] B. Kågström and A. Ruhe, An algorithm for numerical computation of the Jordan normal form of a complex matrix, ACM Trans. Math. Software 6 (1980), 398-419.

[13] T. Kato, Perturbation Theory for Linear Operators, Grundlehren der mathematischen Wissenschaften 132, Springer-Verlag, New York (1980).

[14] P. Lancaster and M. Tismenetsky, The Theory of Matrices, Academic Press, New York, New York (1985).

[15] L. V. Levantovskii, On singularities of the boundary of a stability region, Moscow University Math. Bull. 35:6(1980), 19-22.

[16] L. V. Levantovskii, The boundary of a set of stable matrices, Russian Math. Surveys 35:2(1980), 249-250.

[17] L. V. Levantovskii, Singularities of the boundary of the stability domain, Funct. Anal. Appl. 16 (1982), 34-37.

[18] Matlab User’s Guide, The Mathworks, Sherborn, MA (1989).

[19] Sir I. Newton, Enumeration of Lines of Third Order, Generation of Curves by Shadows, Organic Description of Curves, and Construction of Equations, Translation from Latin by C. R. M. Talbot, Printed by H. G. Bohn, York Street, Covent Garden, London (1860).

[20] M. L. Overton, Large-scale optimization of eigenvalues. SIAM J. Optimization 3 (1992), to appear.

[21] M. L. Overton and R. S. Womersley, On minimizing the spectral radius of a nonsymmetric matrix function: Optimality conditions and duality theory, SIAM J. Matrix Anal. Appl. 9(1988), 473-498.

[22] V. Puiseux, Recherches sur les fonctions algébriques, J. Math. Pures Appl. 51(1850), 465-581.

[23] M. M. Vainberg and V. A. Trenogen, Theory of Branching of Solutions of Non-Linear Equations, Noordhoff International Publishing, Leyden, The Netherlands (1974).

[24] J.H. Wilkinson, The Algebraic Eigenvalue Problem, Clarendon Press, Oxford (1965).