On the Subdifferentiability of Functions of a Matrix Spectrum I: Mathematical Foundations

James V. Burke^{∗} and Michael L. Overton^{†}

Abstract

We consider analytic matrix valued mappings A: C ↦ C^{n×n} and study the variational properties of the spectrum of A(ϵ). Of particular interest are α(ϵ) and ρ(ϵ), respectively the spectral abscissa and the spectral radius of A(ϵ). In this paper, we introduce the mathematical techniques required for this analysis. We begin with polynomials and discuss the bifurcation of the roots of a polynomial having analytic coefficients. It is this bifurcation phenomenon that leads to the nonlipschitzian behavior of the types of functions that we wish to study. Puiseux-Newton series and diagrams are then introduced as a means for analyzing these bifurcations. It is shown how these techniques can be used to describe the tangent cone to certain sets of stable polynomials.

Matrices and polynomials are connected via characteristic polynomials. Further properties of the spectrum of a matrix A^{(0)} are obtained from a block diagonalization of A^{(0)}, where the kth diagonal block is an n_{k} by n_{k} upper triangular matrix, with a constant diagonal whose value is an eigenvalue of A^{(0)} with multiplicity n_{k}. By using results of Arnold on the versal deformation of a matrix, we show how the results concerning polynomials can be translated into into results about matrices. In the case that the block diagonal form is the Jordan form, the conditions which are obtained reduce to conditions on generalized traces of a matrix.

1 Introduction

In this study we consider the variational properties of two related functions of the spectrum of an analytic matrix valued function. Let A be an analytic matrix valued mapping from C to C^{n×n}; thus each element of A(ϵ) is an analytic (holomorphic) function of a single complex parameter ϵ. Define

∗Department of Mathematics, University of Washington. The work of this author was supported in part by National Science Foundation Grant DMS-88-03206

†Computer Science Department, Courant Institute of Mathematical Sciences, New York University. The work of this author was supported in part by National Science Foundation Grant CCR-88-02408.

α(ϵ) = max{Re λ : λ ∈ Σ(ϵ)},  (1)

and

ρ(ϵ) = max{|λ|: λ ∈ Σ(ϵ)},  (2)

where Σ(ϵ) is the spectrum of A(ϵ), i.e. the roots of the characteristic polynomial

P (ϵ, λ) = det(λI − A(ϵ)) = 0.  (3)

The elements of Σ(ϵ) are called the eigenvalues of A(ϵ), the function α(ϵ) is called the spectral abscissa of A(ϵ), and ρ(ϵ) is called the spectral radius. The spectral abscissa and radius are associated with stability properties of the matrix A(ϵ) and are important quantities in many applications. The spectral abscissa is relevant when stability is defined in terms of exponentiation of a matrix, while the spectral radius is relevant if stability is defined in terms of matrix powers. The former typically arises in applications involving differential equations while the latter is relevant in the case of difference equations (which themselves often arise as approximations to differential equations).

Both the functions α(ϵ) and ρ(ϵ) are, in general, nonlipschitzian. The reason for this is the well known phenomenon of bifurcation of roots that occurs when a polynomial with a multiple root is subjected to analytic perturbation. This behavior was well understood by Newton as early as 1680 and lead to his development of so-called Puiseux-Newton series as a means for describing these roots. (The name Puiseux is associated with these series because it was he who, two hundred years later, proved their convergence [10].) The Puiseux-Newton series is a power series in fractional powers of the perturbation parameter ϵ, with the smallest such power being greater than or equal to the inverse of the degree of the polynomial. Thus, in general, α(ϵ) and ρ(ϵ) can only be said to be Hölder continuous of order 1/n. A simple example is:

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

The eigenvalues are the roots of

λ^{n} − ϵ = 0

i.e. the quantity ϵ^{1/n} times the nth roots of unity. Taking ϵ to be real and positive, we thus have α(ϵ) = ρ(ϵ) = ϵ^{1/n}. A more interesting example will be given shortly.

Throughout Part I of this paper, the analytical dependence of the matrix A is restricted to a single complex variable ϵ, which we will often take to be real and positive. In Part II [4] we shall consider dependence of A on several complex variables, and then the results of Part I will be applied to obtain results along curves in the parameter space of Part II.

2 Roots of Polynomials

Matrix theoretic results can always be applied to obtain insight into the behavior of polynomials by stating the results in terms of companion matrices. On the other hand, polynomials are crude measures of the behavior of matrices since their structure is far less rich. Nonetheless, polynomial results provide an important starting point for investigation into matrices.

Because the eigenvalues are the roots of the characteristic polynomial P (ϵ, λ), we begin by addressing the conditions on the coefficients of the characteristic polynomial that are imposed by requiring α(ϵ) or ρ(ϵ) to have certain properties. A key point is to simplify the class of polynomials that one needs to consider. From [2, pp.376-381], there is an open disk in the complex plane containing the origin on which the polynomial P (ϵ, λ) has the unique representation

P (ϵ, λ) = ∏_{λ_{k}∈Σ(0)}((λ_{k} − λ)^{n_{k}} + β_{k1}(ϵ)(λ_{k} − λ)^{n_{k}−1} + ⋯ + β_{kn_{k}}(ϵ))  (4)

where the product index runs over each eigenvalue λ_{k} (with multiplicity n_{k}) in Σ(0), and the β_{kj}(ϵ) are analytic functions vanishing at ϵ = 0. For the moment we shall assume that there is only one such eigenvalue λ_{0}, with multiplicity n_{0}, so that P reduces to

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

with

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

We now develop the key result needed to understand the variational behavior of the spectral abscissa and radius. The result provides information about the first-order behavior of the coefficients of P (ϵ, λ) under an assumption on the rate of change of its roots in a certain direction in the complex plane.

Lemma 1 Consider the polynomial equation (5). Let y_{0} be any nonzero complex number and suppose that all the roots λ(ϵ) of (5) satisfy

Re \overline{y}_{0}(λ(ϵ) − λ_{0}) ≤ δϵ + o(ϵ)  (6)

for some real sequence {ϵ^{ν}} with ϵ^{ν} ↓ 0. Then

Re \overline{y}_{0}β_{1}^{(1)} ≥ −n_{0}δ,  (7)
Re \overline{y}_{0}^{2}β_{2}^{(1)} ≥ 0, Im \overline{y}_{0}^{2}β_{2}^{(1)} = 0,  (8)
β_{j}^{(1)} = 0, j = 3,…, n_{0}.  (9)

An important special case is y_{0} = 1, δ = 0. In this case, Lemma 1 states the following: the tangent cone of the set of stable polynomials, i.e. those with roots in the left half-plane, at the point (λ − λ_{0})^{n} in the space of monic polynomials with analytic coeffients, is contained in the set of polynomials (5) satisfying Re β_{1}^{(1)} ≥ 0, Re β_{2}^{(1)} ≥ 0, Im β_{2}^{(1)} = 0, and β_{j}^{(1)} = 0, j = 3,…, n_{0}. In fact, this set equals the tangent cone [3,8].

The proof of this lemma may be found in [5]. It uses the Puiseux-Newton diagram, a technique devised by Newton [9] for computing the coefficients of Puiseux-Newton series. Although we shall not repeat the proof of Lemma 1 here, we shall give some motivating remarks by means of an example.

Suppose that λ_{0} = 0 and

P (ϵ, λ) = λ^{3} + ϵλ^{2} + (−ϵ − ϵ)^{2}λ + (ϵ^{2} + 2ϵ^{3}).  (10)

Consider the diagram in Figure 1, in which (a) the power of ϵ in the leading nonzero term of the jth coefficient is plotted against j (this includes the point (0, 0)), and (b) the lower boundary of the convex hull of these points, namely a piecewise linear function, is drawn:

[figure: j; X; X; X; X; 1; 2; 3; 2; 1]

Figure 1

Newton used an “Ansatz” argument to show that the slopes of this piecewise linear function are precisely the powers of ϵ in the leading terms of the expansions for λ(ϵ), the roots of P (ϵ, λ) = 0. For this example, the roots have the form

λ(ϵ) = ±ϵ^{\frac{1}{2}} + ⋯,  (11)

and

λ(ϵ) = ϵ + ⋯,  (12)

as can be verified by substitution into (10) and observing cancellation. Note that there are two roots (11) whose leading power of ϵ is 1/2, reflecting the fact that the line segment in Figure 1 with slope 1/2 runs from j = 0 to j = 2. The coefficients of these ϵ^{\frac{1}{2}} terms are, in this case, the two square roots of −β_{2}^{(1)} = 1; if β_{2}(ϵ) is changed to ϵ − ϵ^{2}, we obtain leading terms ±iϵ^{\frac{1}{2}} in (11), where i = \sqrt{−1}. Similarly, there is only one root (12) with leading power of ϵ equal to 1, reflecting the fact that the line segment with slope 1 runs from j = 2 to j = 3, and the coefficient of this leading term is −β_{3}^{(2)}/β_{2}^{(1)} = 1. (That these are the relevant coefficients follows from the fact that the line segment with slope 1 interpolates the points (2,1) and (3,2).)

With this example understood, the result of Lemma 1 can now be explained. Let y_{0} have any nonzero value, and let L be the line in the complex plane defined by \{z : Re \overline{y}_{0}(z −λ_{0}) = 0}. In order for (6) to hold, all slopes in the associated Puiseux-Newton diagram must be ≥ 1/2 since, for example, a slope with value 1/3 corresponds to three roots with leading term ϵ^{\frac{1}{3}} and coefficients equal to the three cube roots of some complex number, which means that at least one of the roots must lie on one side of L and at least one of the roots must lie on the other side. This, then, explains (9). In the case of a slope with value 1/2, the only possible way (6) can hold is if both roots lie on L; this amounts to a condition on β_{2}^{(1)} which is given by (8). Finally, (7) follows directly from (6) since −β_{1}(ϵ) is the sum of the roots λ(ϵ) (shifted by λ_{0}).

The result of Lemma 1 extends immediately to general polynomials of the form (4). However, we shall not need to state this explicitly; instead we go on to consider the original matrix problem.

3 Eigenvalues of Matrices

A matrix can be reduced via similarity transformations to a variety of canonical forms. A finite number of elementary unitary transformations is sufficient to reduce a matrix to Hessenberg form, where all subdiagonals except the first are reduced to zero; this can be further reduced to upper triangular or Schur form by a general unitary transformation. The spectrum of a matrix appears on the diagonal of its Schur form, but other information, such as invariant subspace information, is not apparent. In order to further reduce the matrix, i.e. to introduce zeros in the upper triangle as well as the lower, nonunitary transformations are generally required; such transformations are potentially quite ill-conditioned, i.e. the norm of the transformation times the norm of its inverse could be large. The ultimate canonical form is the Jordan form, where zeros are introduced everywhere except on the diagonal and some parts of the first superdiagonal. However, the Jordan form is a discontinuous function over the space of matrices. See [7] for an excellent general discussion.

The canonical form that we shall need is block diagonal form ([7, Section 7.1.3]). This generally requires nonunitary transformations but is not as difficult to compute as the Jordan form (which is a special case). Let A(ϵ) be the analytic matrix valued function of Section 1, and assume that

A(0) = A^{(0)} = P DP^{−1}, D = Diag(D_{1},…, D_{m}),  (13)

where the kth diagonal block N_{k} is upper triangular with constant diagonal, i.e.

N_{k} = λ_{k}I + N_{k},

with N_{k} strictly upper triangular and hence nilpotent. Here m is the number of distinct eigenvalues of A, so λ_{k} is an eigenvalue of A(0) with (algebraic) multiplicity n_{k}, the order of N_{k}. The geometric multiplicity of λ_{k} is the nullity of N_{k}; this is easily seen to be the number of independent eigenvectors associated with λ_{k}. If N_{k} = 0, the eigenvalue λ_{k} is said to be semisimple (or nondefective); if N_{k} has rank n_{k} − 1, λ_{k} is said to be nonderogatory.

Now define G(ϵ) = P^{−1}A(ϵ)P, so that G(0) = D. A result of V.I. Arnold [1] states that G(ϵ) has the following versal deformation:

G(ϵ) = Y (ϵ)H(ϵ)Y (ϵ)^{−1}  (14)

where Y and H are both analytic, Y (0) = I, and H(ϵ) commutes with D for all ϵ in a neighborhood of 0. Now since D is block diagonal, H(ϵ) must also be block diagonal with the same block sizes n_{1},…, n_{m} (one way to prove this is that D is similar, via a block diagonal similarity transformation, to its Jordan form, and the matrices commuting with a Jordan form have a special block diagonal structure [1,6]). Now A(ϵ), G(ϵ) and H(ϵ) are all similar to each other, so the eigenvalues of A(ϵ) are given by the roots of the characteristic polynomial of H(ϵ), which is the product of characteristic polynomials of the diagonal blocks H_{k}(ϵ). The following lemma, brought to our attention by J. Sylvester[11], shows how to compute the derivatives of the coefficients of these characteristic polynomials at ϵ = 0:

Lemma 2

\frac{d}{dϵ} det(λI − H_{k}(ϵ))|_{ϵ=0}  =  −tr (H_{k}^{(1)})(λ − λ_{k})^{n_{k}−1} − tr (N_{k}H_{k}^{(1)})(λ − λ_{k})^{n_{k}−2}
− ⋯ − tr (N_{k}^{n_{k}−2}H_{k}^{(1)})(λ − λ_{k}) − tr (N_{k}^{n_{k}−1}H_{k}^{(1)}),

where H_{k}^{(1)} = H′_{k}(0).

Proof First note that H_{k}(0) is the kth diagonal block of G(0), i.e. D_{k}, so we may write H_{k}(ϵ) = D_{k} + ϵH_{k}^{(1)} + o(ϵ). Let μ = λ − λ_{k}, so

λI − H_{k}(ϵ) = μI − N_{k} − ϵH_{k}^{(1)} + o(ϵ).

Then

\frac{d}{dϵ} det(λI − H_{k}(ϵ))|_{ϵ=0}  =  det(μI − N_{k})\frac{d}{dϵ} det(I − ϵμ^{−1}(I − μ^{−1}N_{k})^{−1}H_{k}^{(1)})
=  −μ^{n_{k}}tr (μ^{−1}(I − μ^{−1}N_{k})^{−1}H_{k}^{(1)})
= −μ^{n_{k}−1} tr (I + μ^{−1}N_{k} + (μ^{−1}N_{k})^{2} + ⋯ + (μ^{−1}N_{k})^{n_{k}−1})H_{k}^{(1)}.

We are now in a position to prove the main result of this section.

Lemma 3 Let A(ϵ) = A^{(0)} +ϵA^{(1)} +⋯ be an analytic matrix function with A^{(0)} having the block diagonal decomposition (13), and let B_{1},…, B_{m} be the corresponding diagonal blocks of P^{−1}A^{(1)}P (this matrix is not block diagonal in general). Then the eigenvalues λ(ϵ) of A(ϵ) corresponding to the eigenvalue λ_{k} of A^{(0)} are the roots of the polynomial

(λ − λ_{k})^{n_{k}} + β_{k1}(ϵ)(λ − λ_{k})^{n_{k}−1} + ⋯ + β_{kn_{k}}(ϵ) = 0

where β_{kj} are analytic functions with

β′_{kj}(0) = −tr N_{k}^{j−1}B_{k}.

Proof Using the versal deformation (14), we have

P^{−1}A′(0)P = G′(0) = H′(0) + Y′(0)H(0) − H(0)Y′(0),

i.e.

P^{−1}A^{(1)}P = H^{(1)} + Y^{(1)}D − DY^{(1)}.

Therefore

tr (N_{k}^{j−1}B_{k}) = tr (N_{k}^{j−1}H_{k}^{(1)}) + tr (N_{k}^{j−1}(Y_{k}^{(1)}D_{k} − D_{k}Y_{k}^{(1)})).

where Y_{k}^{(1)} is the kth diagonal block of Y′(0). Now note that N_{k}^{j−1} commutes with D_{k}, so the last term is zero. (This follows since the trace of E(ZD − DZ) is zero if and only if E commutes with D; this easily verified fact is a basic tool in the derivation of the Arnold versal deformation.) The proof is completed by the application of Lemma 2.

Remark. Suppose that the block decomposition is the Jordan form, i.e. all superdiagonals of N_{k} are zero except the first, which consists of zeros and ones. Suppose further that all eigenvalues are nonderogatory, i.e. no first superdiagonal contains a zero. Then the quantity tr N_{k}^{j−1}B_{k} reduces to tr ^{(j)}B_{k}, the jth generalized trace of B_{k}, which is defined to be the sum of the elements on the (j − 1)th subdiagonal of B_{k}. In the derogatory case, we obtain a sum of such generalized traces.

4 The Spectral Abscissa

We now obtain the following result:

Lemma 4 Let A(ϵ) = A^{(0)} +ϵA^{(1)} +⋯ be an analytic matrix function with A^{(0)} having the block diagonal decomposition (13), and let B_{1},…, B_{m} be the corresponding diagonal blocks of P^{−1}A^{(1)}P. Recall the definition of the spectral abscissa α(ϵ) in (1), and define

A = {k : Re λ_{k} = α(0)}.

Suppose that

α(ϵ) − α(0) ≤ δϵ + o(ϵ)  (15)

for some real sequence {ϵ^{ν}} with ϵ^{ν} ↓ 0. Then, for each k ∈ A,

Re tr B_{k} ≤ n_{k}δ,
Re tr (N_{k}B_{k}) ≤ 0, Im tr (N_{k}B_{k}) = 0,
tr (N_{k}^{j−1}B_{k}) = 0, j = 3,…, n_{k}.

Proof Clearly, (15) is equivalent to requiring, for each k ∈ A, that

Re (λ(ϵ) − λ_{k}) ≤ δϵ + o(ϵ)

for all λ(ϵ) which are eigenvalues of A(ϵ) corresponding to λ_{k}. The proof therefore follows from applying Lemma 1 (with y_{0} = 1) and Lemma 3.

5 The Spectral Radius

Similarly, we obtain

Lemma 5 Using the same assumptions as in the previous lemma, recall the definition of the spectral radius ρ(ϵ) in (2), and define

R = {k :|λ_{k}|= ρ(0)}.

Suppose that

ρ(ϵ) − ρ(0) ≤ δϵ + o(ϵ)  (16)

for some real sequence {ϵ^{ν}} with ϵ^{ν} ↓ 0. Then, for each k ∈ R,

Re (\overline{λ}_{k}tr B_{k})+ |tr N_{k}B_{k}|≤ n_{k}δ |λ_{k}|,  (17)
Re (\overline{λ}_{k}^{2}tr (N_{k}B_{k})) ≤ 0, Im (\overline{λ}_{k}^{2}tr (N_{k}B_{k})) = 0,  (18)
tr (N_{k}^{j−1}B_{k}) = 0, j = 3,…, n_{k}.

Partial Proof Equation (16) is equivalent to requiring, for each k ∈ R, that

|λ(ϵ)| − |λ_{k}|≤ δϵ + o(ϵ)

for all λ(ϵ) which are eigenvalues of A(ϵ) corresponding to λ_{k}, i.e.

|λ(ϵ)|^{2} − |λ_{k}|^{2}≤ δϵ(|λ(ϵ)| + |λ_{k}|) + o(ϵ)

or equivalently

Re (\overline{λ}_{k}(λ(ϵ) − λ_{k})) + \frac{1}{2} |λ(ϵ) − λ_{k}|^{2}≤ \frac{1}{2}δϵ(|λ(ϵ)| + |λ_{k}|) + o(ϵ).

Dropping a positive term from the left-hand side, and using the Hölder continuity of λ(ϵ), we have

Re (\overline{λ}_{k}(λ(ϵ) − λ_{k})) ≤ δϵ |λ_{k}| +o(ϵ).

By applying Lemma 1, with y_{0} = λ_{k}, and Lemma 3, we obtain almost the result we need, but without the second term on the left-hand side of (17). Obtaining (17) requires a more careful argument, which is given in [5].

References

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

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

[3] J. V. Burke and M. L. Overton, Stable perturbations of nonsymmetric matrices, Computer Science Dept Report 520, New York University, Oct 1990.

[4] J. V. Burke and M. L. Overton, On the subdifferentiability of functions of a matrix spectrum II: Subdifferential forumulas, this volume.

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

[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] G.H. Golub and C. Van Loan, Matrix Computations, Johns Hopkins University Press, 2nd edition (1989).

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

[9] 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).

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

[11] J. Sylvester, University of Washington. Private communication.