Äîêóìåíò âçÿò èç êýøà ïîèñêîâîé ìàøèíû. Àäðåñ îðèãèíàëüíîãî äîêóìåíòà : http://mailybaev.imec.msu.ru/papers/Mailybaev2005.pdf Äàòà èçìåíåíèÿ: Tue Jun 14 12:10:04 2005 Äàòà èíäåêñèðîâàíèÿ: Mon Oct 1 19:37:51 2012 Êîäèðîâêà:

Computation of multiple eigenvalues and generalized eigenvectors for matrices dependent on parameters
arXiv:math-ph/0502010 v1 2 Feb 2005
Alexei A. Mailybaev


Abstract The pap er develops Newton's method of finding multiple eigenvalues with one Jordan block and corresp onding generalized eigenvectors for matrices dep endent on parameters. It computes the nearest value of a parameter vector with a matrix having a multiple eigenvalue of given multiplicity. The method also works in the whole matrix space (in the absence of parameters). The approach is based on the versal deformation theory for matrices. Numerical examples are given.

Keywords: multiparameter matrix family, multiple eigenvalue, generalized eigenvector, Jordan blo ck, versal deformation, Schur decomposition

1

Introduction

Transformation of a square nonsymmetric (non-Hermitian) matrix A to the Jordan canonical form is the classical sub ject that finds various applications in pure and applied mathematics and natural sciences. It is well known that a generic matrix has only simple eigenvalues and its Jordan canonical form is a diagonal matrix. Nevertheless, multiple eigenvalues typically appear in matrix families, and one Jordan blo ck is the most typical Jordan structure of a multiple eigenvalue [3, 4]. Many interesting and important phenomena asso ciated with qualitative changes in the dynamics of mechanical systems [20, 29, 30, 36], stability optimization [6, 21, 25], and bifurcations of eigenvalues under matrix perturbations [32, 35, 34, 38] are related to multiple eigenvalues. Recently, multiple eigenvalues with one Jordan blo ck became of great interest in physics, including quantum mechanics and nuclear physics [2, 17, 24], optics [5], and electrical engineering [8]. In most applications, multiple eigenvalues appear through the intro duction of parameters. In the presence of multiple eigenvalues, the numerical problem of computation of the Jordan canonical form is unstable, since the degenerate structure can be destroyed by arbitrarily small perturbations (caused, for example, by round-off errors). Hence, instead of analyzing a single matrix, we should consider this problem in some neighborho o d in matrix or parameter space. Such formulation leads to the important problem left open
Institute of Mechanics, Lomonosov Moscow State University, Michurinsky pr. 1, 119192 Moscow, Russia. E-mail: mailybaev@imec.msu.ru, Tel.: (7095) 939 2039, Fax: (7095) 939 0165.


1


by Wilkinson [40, 41]: to find the distance of a given matrix to the nearest degenerate matrix. We study the problem of finding multiple eigenvalues for matrices dependent on several parameters. This implies that matrix perturbations are restricted to a specific submanifold in matrix space. Such restriction is the main difficulty and difference of this problem from the classical analysis in matrix spaces. Existing approaches for finding matrices with multiple eigenvalues [7, 9, 11, 12, 16, 18, 19, 23, 26, 33, 40, 41] assume arbitrary perturbations of a matrix and, hence, they do not work for multiparameter matrix families. We also mention the topological metho d for the lo calization of double eigenvalues in twoparameter matrix families [22]. In this paper, we develop Newton's metho d for finding multiple eigenvalues with one Jordan blo ck and corresponding generalized eigenvectors in multiparameter matrix families. The presented metho d solves numerically the Wilkinson problem of finding the nearest matrix with a multiple eigenvalue (both in multiparameter and matrix space formulations). The implementation of the metho d in MATLAB co de is available, see [31]. The metho d is based on the versal deformation theory for matrices. In spirit, our approach is similar to [13], where matrices with multiple eigenvalues where found by path-following in matrix space (multiparameter case was not considered). The paper is organized as follows. In Section 2, we intro duce concepts of singularity theory and describe a general idea of the paper. Section 3 provides expressions for values and derivatives of versal deformation functions, which are used in Newton's metho d in Section 4. Section 5 contains examples. In Section 6 we discuss convergence and accuracy of the metho d. Section 7 analyzes the relation of multiple eigenvalues with sensitivities of simple eigenvalues of perturbed matrices. In Conclusion, we summarize this contribution and discuss possible extensions of the metho d. Pro ofs are collected in the Appendix.

2

Multiple eigenvalues with one Jordan block in multiparameter matrix families

Let us consider an m â m complex non-Hermitian matrix A, which is an analytical function of a vector of complex parameters p = (p1 , . . . , pn ). Similarly, one can consider real or complex matrices smo othly dependent on real parameters, and we will comment the difference among these cases where appropriate. Our goal is to find the values of parameter vector p at which the matrix A(p) has an eigenvalue of algebraic multiplicity d with one d â d Jordan blo ck (geometric multiplicity 1). Such and eigenvalue is called nonderogatory. There is a Jordan chain of generalized vectors u1 , . . . , ud (the eigenvector and asso ciated vectors) corresponding to and determined by the equations Au1 = u1 , Au2 = u2 + u1 , . . . Aud = ud + ud-1 . 2

(2.1)


Figure 1: Geometry of the bifurcation diagram. These vectors form an m â d matrix U = [u1 , . . . , ud ] satisfying the equation 1 . .. AU = UJ , J = , .. . 1

(2.2)

where J is the Jordan blo ck of size d. Recall that the Jordan chains taken for all the eigenvalues and Jordan blo cks determine the transformation of the matrix A to the Jordan canonical form [14]. In singularity theory [4], parameter space is divided into a set of strata (smo oth submanifolds of different dimensions), which correspond to different Jordan structures of the matrix A. Consider, for example the matrix family 0100 p 0 1 0 (2.3) A(p) = 1 p2 0 0 1 , p = (p1 , p2 , p3 ). p3 0 0 0 The bifurcation diagram in parameter space is shown in Figure 1 (for simplicity, we consider only real values of parameters). There are four degenerate strata: 2 (surfaces), 3 and 2 2 (curves), and 4 (a point). The surface 2 , curve 3 , and point 4 corre12 spond, respectively, to the matrices with double, triple, and quadruple eigenvalues with one Jordan blo ck. The curve 2 2 is the transversal self-intersection of the stratum 2 12 corresponding to the matrices having two different double eigenvalues. This bifurcation diagram represents the well-known "swallow tail" singularity [4]. We study the set of parameter vectors, denoted by d , corresponding to matrices having multiple eigenvalues with one Jordan blo ck of size d. The set d is a smo oth surface in parameter space having co dimension d - 1 [3, 4]. Thus, the problem of finding multiple eigenvalues in a matrix family is equivalent to finding the surface d or its particular point. Since the surface d is smo oth, we can find it numerically by using Newton's metho d. This requires describing the surface d as a solution of d - 1 equations qi (p) = 0, 3 i = 2, . . . , d , (2.4)


for independent smo oth functions qi (p). (In these notations, we the multiple eigenvalue = q1 (p).) Finding the functions qi (p) is the clue to the problem solution. In this paper, we define the functions qi (p) in the following deformation theory [3, 4], in the neighborho o d of d , the matrix q1 (p) 1 q (p) q1 (p) A(p)U(p) = U(p)B(p), B(p) = 2 . . . qd (p)

keep the first function for and their first derivatives way. According to versal A(p) satisfies the relation .. . .. . 1 q1 (p) ,

(2.5)

where U(p) is an m â d analytic matrix family, and q1 (p), . . . , qd (p) are analytic functions (blank places in the matrix are zeros). The functions q1 (p), . . . , qd (p) are uniquely determined by the matrix family A(p). By using (2.5), it is straightforward to see that the surface d is defined by equations (2.4). If (2.4) are satisfied, the matrix B(p) is the d â d Jordan blo ck. Hence, at p d , the multiple eigenvalue is = q1 (p) and the columns of U(p) are the generalized eigenvectors satisfying equations (2.1). The metho d of finding the functions qi (p) and U(p) and their derivatives at the point p d has been developed in [27, 28]. In Newton's metho d for solving (2.4), we need the values and derivatives of the functions qi (p) at an arbitrary point p d . /

3

Linearization of versal deformation functions

Let p0 be a given parameter vector determining a matrix A0 = A(p0 ). Since multiple eigenvalues are nongeneric, we typically deal with a diagonalizable matrix A0 . Let 1 , . . . , m be eigenvalues of the matrix A0 . We sort these eigenvalues so that the first d of them, 1 , . . . , d , coalesce as the parameter vector is transferred continuously to the surface d . The eigenvalues that form a multiple eigenvalue are usually known from the context of a particular problem. Otherwise, one can test different sets of d eigenvalues. Let us cho ose m â d matrices X and Y such that A0 X = XS, Y A0 = SY , Y X = I,
d

(3.1)

where S is the d â d matrix whose eigenvalues are 1 , . . . , conjugate transpose. The first two equalities in (3.1) imply X span the right invariant subspace of A0 corresponding of Y span the left invariant subspace. The third equality The matrix S can be expressed as S = Y A0 X,

; the star denotes the complex that the columns of the matrix to 1 , . . . , d , and the columns is the normalization condition. (3.2)

which means that S is the restriction of the matrix operator A0 to the invariant subspace given by the columns of X. The constructive way of cho osing the matrices S, X, and Y will be described in the next section. 4


The following theorem provides the values and derivatives of the functions q1 (p), . . . , qd (p) in the versal deformation (2.5) at the point p0 . Theorem 3.1 Let S, Y , and X be the matrices satisfying equations (3.1). Then q1 (p0 ) = trace S/d, (3.3)

and the values of q2 (p0 ), . . . , qd (p0 ) are found as the characteristic polynomial coefficients of the traceless matrix S - q1 (p0 )I: z d - q2 (p0 )z
d-2

- · · · - qd-1 (p0 )z - qd (p0 ) = det (z + q1 (p0 ))I - S ,

(3.4)

where I is the d â d identity matrix. The first derivatives of the functions qi (p) at p0 are determined by the recurrent formulae A q1 = trace Y X /d , pj pj qi = trace (S - q1 (p0 )I) pj i = 2, . . . , d , j = 1, . . . , n,
i-1

A Y X - trace(C pj


i-1

q1 ) - pj

i-1

trace(C
k =2

i-1

Ek 1 )

qk , pj

(3.5) where the derivatives are evaluated at p0 ; C = B(p0 ) - q1 (p0 )I is the companion matrix 0 1 d . q2 (p0 ) 0 . . C= qi (p0 )Ei1 , (3.6) = J0 + . .. . . 1 . i=2 qd (p0 ) 0 and Ei1 is the matrix having the unit (i, 1)th element and zeros in other places. The pro of of this theorem is given in the Appendix. When the matrix A is arbitrary (not restricted to a multiparameter each entry of the matrix can be considered is an independent paramet matrix A can be used instead of the parameter vector: p - A. The with respect to its (i, j )th entry is Eij . Thus, the formulae of Theorem 3.1 matrix family), er. Hence, the derivative of A can be applied.

Corollary 3.1 Let S, Y , and X be the matrices satisfying equations (3.1). Then the values of q1 (A0 ), . . . , qd (A0 ) are given by formulae (3.3) and (3.4) with p0 substituted by A0 . Derivatives of the functions q1 (A), . . . , qd (A) with respect to components of the matrix A taken at A0 are q1 = (XY /d)T , A qi = X(S - q1 (p0 )I) A i = 2, . . . , d .
i-1

Y

T

- trace(C

i-1

)

q1 - A

i-1

trace(C
k =2

i-1

Ek 1 )

qk , A

(3.7)

5


Here T is the transpose operator, and

··· qi .. = . A qi ··· am1

qi a11 . . .

qi a1m . . . qi amm

(3.8)

is the m â m matrix of derivatives of qi (A) with respect to components of the matrix A taken at A0 . At p0 d , we can find the multiple eigenvalue and the corresponding Jordan chain of generalized eigenvectors u1 , . . . , ud . This problem reduces to the transformation of the matrix S to the prescribed Jordan form (one Jordan blo ck). A possible way of solving this problem is presented in the following theorem (see the Appendix for the pro of.). Theorem 3.2 At the point p0 d , the multiple eigenvalue is given by the expression = trace S/d. The general form of the Jordan chain of generalized eigenvectors u1 , . . . , ud is u1 = X(S - I)
d-1

(3.9)

k,

...,

ud-1 = X(S - I)k, ud = Xk,

(3.10)

where k Cd is an arbitrary vector such that the eigenvector u1 is nonzero. Choosing a ^ particular unit-norm eigenvector u1 , e.g., by taking the scaled biggest norm column of the d-1 matrix X(S - I) , one can fix the vector k by the orthonormality conditions ^1 u ui = 1, i = 1; 0 , i = 2, . . . , d . (3.11)

The accuracy of the multiple eigenvalue and generalized eigenvectors determined by formulae (3.9) and (3.10) has the same order as the accuracy of the point p0 in the surface d .

4

Newton's method

There are several ways to find the matrices S, X, and Y . The simplest way is to use the diagonalization of A0 . Then S = diag(1 , . . . , d ) is the diagonal matrix, and the columns of X and Y are the right and left eigenvectors corresponding to 1 , . . . , d . This way will be discussed in Section 5. If the parameter vector p0 is close to the surface d , the diagonalization of the matrix A0 is ill-conditioned. Instead of the diagonalization, one can use the numerically stable Schur decomposition S = X A0 X, where S is an upper-triangular matrix called the Schur canonical form, and X = (X )-1 is a unitary matrix [15]. The diagonal elements 6


s11 , . . . , smm of S are the eigenvalues of A0 . We can cho ose the Schur form so that the first d diagonal elements are the eigenvalues 1 , . . . , d . Performing the blo ck-diagonalization of the Schur form S [15, §7.6], we obtain the blo ck-diagonal matrix S0 0 S = [Y , Y ] A0 [X, X ], (4.1)

where S is a d â d upper-triangular matrix with the diagonal (1 , . . . , d ); [X, X ] and [Y , Y ] = [X, X]-1 are nonsingular m â m matrices (not necessarily unitary). These operations with a Schur canonical form are standard and included in many numerical linear algebra packages, for example, LAPACK [1]. They are numerically stable if the eigenvalues 1 , . . . , d are separated from the remaining part of the spectrum. As a result, we obtain the matrices S, X, and Y satisfying equations (3.1). When the matrices S, X, and Y are determined, Theorem 3.1 provides the necessary information for using Newton's metho d for determining the stratum d . Indeed, having the parameter vector p0 = (p0 , . . . , p0 ) as the initial guess, we linearize equations (2.4) of 1 n the surface d as n qi qi (p0 ) + (p j - p 0 ) = 0 , i = 2 , . . . , d , (4.2) j pj j =1 where the values of qi (p0 ) and the derivatives qi / pj at p0 are provided by Theorem 3.1. In the generic case, the linear part in (4.2) is given by the maximal rank matrix [ qi / pj ]. System (4.2) has the single solution if the number of parameters n = d - 1 (the set d is an isolated point). If n > d - 1, one can take the least squares solution or any other solution depending on which point of the surface d one would like to find. If n < d, the multiple eigenvalue still can exist in matrices with symmetries (e.g., Hamiltonian or reversible matrices [35]); then the least squares fit solution of (4.2) is a go o d choice. In Newton's metho d, the obtained vector of parameters p = (p1 , . . . , pn ) is used in the next iteration. In each iteration, we should cho ose d eigenvalues of the matrix A0 . These are the d eigenvalues nearest to the approximate multiple eigenvalue
n

= q1 (p0 ) +
j =1

q1 (p j - p 0 ) j pj

(4.3)

calculated at the previous step. If the iteration pro cedure converges, we obtain a point p d . Then the multiple eigenvalue and corresponding generalized eigenvectors are found by Theorem 3.2. Note that, at the point p d , system (4.2) determines the tangent plane to the surface d in parameter space. The pseudo-co de of the described iteration pro cedure is presented in Table 1. Depending on a particular application, the line 3 in this pseudo-co de can be implemented in different ways, e.g., as the least squares solution or as the solution nearest to the input parameter vector p0 . The implementation of this metho d in MATLAB co de is available, see [31]. In case of complex matrices dependent on real parameters, the same formulae can be used. In this case, system (4.2) represents 2(d - 1) independent equations (each equality 7


INPUT: matrix family A(p), initial parameter vector p0 , and eigenvalues 1 , . . . , 1 2 3 4 5 6 7 8 9 : : : : : : : : :

d

Schur decomposition and blo ck-diagonalization (4.1) of the matrix A0 = A(p0 ); evaluate qi (p0 ) and qi / pj by formulae (3.3)­(3.5); find pnew by solving system (4.2) (e.g. the least squares solution); IF pnew - p0 > desired accuracy evaluate approximate multiple eigenvalue app by (4.3); n cho ose d eigenvalues 1 ew , . . . , new of Anew = A(pnew ) nearest to app ; d n perform a new iteration with p0 = pnew and i = i ew , i = 1, . . . , d (GOTO 1); ELSE (IF pnew - p0 desired accuracy) find multiple eigenvalue and generalized eigenvectors by formulae (3.9)­(3.11);

OUTPUT: parameter vector p d , multiple eigenvalue and Jordan chain of generalized eigenvectors u1 , . . . , ud Table 1: Pseudo-co de of Newton's metho d for finding multiple eigenvalues in multiparameter matrix families. determines two equations for real and imaginary parts). This agrees with the fact that the co dimension of d in the space of real parameters is 2(d - 1) [4]. Finally, consider real matrices smo othly dependent on real parameters. For complex multiple eigenvalues, the system (4.2) contains 2(d - 1) independent real equations (co dimension of d is 2(d - 1)). Remark that imaginary parts of the eigenvalues 1 , . . . , d should have the same sign. For real multiple eigenvalues, qi (p0 ) and qi / pj are real (the real Schur decomposition must be used). Hence, (4.2) contains d - 1 real equations (co dimension of d is d - 1). In this case, the eigenvalues 1 , . . . , d are real or appear in complex conjugate pairs. In some applications, like stability theory [35], we are interested in specific multiple eigenvalues, e.g., zero and purely imaginary eigenvalues. In this case equation (4.3) should be included in the linear system of Newton's approximation (4.2). For arbitrary matrices A = [aij ] (without parameters), similar Newton's iteration pro cedure is based on Corollary 3.1. The linearized equations (4.2) are substituted by
m

qi (p0 ) +
j , k =1

qi (aj k - a0k ) = 0, j aj k

i = 2, . . . , d ,

(4.4)

where A0 = [a0k ] is the matrix obtained at the previous step or the initial input matrix. j The first-order approximation of the multiple eigenvalue (4.3) takes the form
m

= q1 (p0 ) +
j , k =1

q1 (aj k - a0k ). j aj k

(4.5)

8


5

Examples

All calculations in the following examples were performed using MATLAB co de [31]. For the sake of brevity, we will show only first several digits of the computation results.

5.1

Example 1

Let us consider the two-parameter family of real matrices 130 A(p) = p1 1 p2 , p = (p1 , p2 ). 231

(5.1)

Bifurcation diagram for this matrix family is found analytically by studying the discriminant of the characteristic polynomial. There are two smo oth curves 2 and a point 3 at the origin (the cusp singularity), see Figure 2. Let us consider the point p0 = (-0.03, 8.99), where the matrix A0 has the eigenvalues 1,2 = -1.995 ± i0.183 and 3 = 6.990. In order to detect a double real eigenvalue, we cho ose the pair of complex conjugate eigenvalues 1 , 2 . By ordering diagonal blo cks in the real Schur form of A0 and blo ck-diagonalizing, we find the matrices S, X, and Y satisfying (3.1) in the form 0.688 -0.676 0.729 -0.574 -1.995 -5.083 , X = -0.688 -0.491 , Y = -0.604 -0.286 . S= 0.007 -1.995 0.231 0.550 0.357 0.858 Applying the formulae of Theorem 3.1, we find q1 (p0 ) q2 (p0 ) = -1.995 -0.033 , q1 / p q2 / p
1 1

q1 / p q2 / p

2 2

=

-0.111 -0.148 1.001 0.333

.

(5.2)

The linearized system (4.2) represents one real scalar equation. We find the nearest parameter vector p 2 (the least squares solution) as p = p0 - q2 (p0 ) ( q2 / p1 , q2 / p2 ) = (-0.00001, 8.99999). ( q2 / p1 )2 + ( q2 / p2 )2 (5.3)

Now let us take different points p0 in the neighborho o d of the curve 2 and calculate one-step Newton's approximations of the nearest points p 2 . In this case we cho ose 1 , 2 as a pair of complex conjugate eigenvalues of A0 = A(p0 ). If all eigenvalues of A0 9

After five iterations of Newton's metho d, we find the Then Theorem 3.2 gives the multiple eigenvalue and 10-15 : 3 1 = -2, [u1 , u2 ] = -3 19 1

exact nearest point p = (0, 9) 2 . the Jordan chain with the accuracy -1 + 30/19 2 - 3 0/1 9 . (5.4) -1 + 10/19


Figure 2: One-step approximations of the nearest points with double eigenvalues. are real, we test all different pairs of eigenvalues, and take the pair providing the nearest point p 2 . The result is shown in Figure 2, where each arrow connects the initial point p0 with the one-step Newton's approximation p. For one point p0 we performed two iterations, taking the point p as a new initial point p0 = p. The convergence of this iteration series is shown in the enlarged part of parameter space (inside the circle in Figure 2). The results confirm Newton's metho d rate of convergence.

5.2

Example 2
0

Let us consider the real matrix A 0 0 A1 = 0

= A1 + E, where 10 342 0 , E = 8 3 6 , 00 496

(5.5)

and = 2.2e-15, = 1.5e-9. This matrix was used in [10] for testing the GUPTRI [18, 19] algorithm. It turned out that this algorithm detects a matrix A 3 (with a nonderogatory triple eigenvalue) at the distance O (10-6) from A0 , while the distance from A0 to 3 is less than E F = 3.62e-14 since A1 3 . This is explained by the observation that the algorithm finds matrix perturbations along a specific set of directions, and these directions are almost tangent to the stratum 3 in the case under consideration [10]. Our metho d determines lo cally the whole stratum 3 in matrix space and, hence, it should work correctly in this case. Since the triple eigenvalue is formed by all eigenvalues of A0 , we can use S = A0 and X = Y = I in the formulae of Corollary 3.1. As a result,

10


we find the least squares solution of system (4.4) in t 0 A = 1.0e-14 -1.760 -0.880 Approximations of the multiple eigenvalue and cor evaluated by Theorem 3.2 for the matrix A are 1.000 0.000 = 8.800e-15, [u1 , u2 , u3 ] = 0.000

he form A = A0 + A, where 00 0 0 . (5.6) 00

responding generalized eigenvectors -0.000 -0.000 1.000 -0.000 . 0.000 6.667e+8

(5.7)

We detected the matrix A at the distance A F = 1.97e-14, which is smaller than the initial perturbation E F = 3.62e-14 ( A F denotes the Frobenius matrix norm). The matrix U = [u1 , u2 , u3 ] satisfies the Jordan chain equation (2.2) with the very high accuracy AU - UJ F / U F = 9.6e-23. The normal complementary subspace N of the tangent space to 3 at A1 has the form [4] 000 N = x 0 0 | x, y R . (5.8) y x 0 It is easy to see that the matrix A in (5.6) is equal to the pro jection of -E to the normal subspace N . This confirms that the obtained matrix A 3 is the nearest to A0 .

5.3

Example 3

Let us consider the 12 â 12 Frank matrix A0 = [a0j ] with the elements i a0j = i n + 1 - max(i, j ), j i - 1, 0, j < i - 1. (5.9)

The Frank matrix has six small positive eigenvalues which are ill-conditioned and form nonderogatory multiple eigenvalues of multiplicities d = 2, . . . , 6 for small perturbations of the matrix. The results obtained by Newton's metho d with the use of Corollary 3.1 are presented in Table 2. An eigenvalue of multiplicity d of the nearest matrix A d is formed by d smallest eigenvalues of A0 . The second column of Table 2 gives the distance dist(A0 , d ) = A - A0 F , where the matrix A is computed after one step of Newton's pro cedure. The third column provides exact distances computed by Newton's metho d, which requires 4­5 iterations to find the distance with the accuracy O (10-15 ). At each iteration, we find the solution A of system (4.4), which is the nearest to the matrix (5.9). The multiple eigenvalues and corresponding generalized eigenvectors are found at the last iteration by Theorem 3.2. The accuracy estimated as AU - UJ F / U F varies between 10-10 and 10-13 . The matrices of generalized eigenvectors U have small condition numbers, which are given in the fourth column of Table 2. For comparison, the fifth and 11


d 2 3 4 5 6

dist(A0 , d ) 1-step approximation 1.619e-10 1.956e-8 1.647e-6 9.299e-5 3.150e-3

dist(A0 , d ) exact 1.850e-10 2.267e-8 1.861e-6 1.020e-4 3.400e-3

cond U 1.125 1.746 4.353 14.14 56.02

A-A

0F

[13]

A-A

0F

[18]

3.682e-10 3.833e-8 3.900e-6 4.280e-4 7.338e-2

6e-3

Table 2: Distances to the multiple eigenvalue strata d for the Frank matrix. sixth columns give upper bounds for the distance to the nearest matrix A d found in [13, 18]. We emphasize that this is the first numerical metho d that is able to find exact distance to a nonderogatory stratum d . Metho ds available in the literature cannot solve this problem neither in matrix space nor for multiparameter matrix families.

6

Convergence and accuracy

In the proposed approach, the standard Schur decomposition and blo ck-diagonalization (4.1) of a matrix are required at each iteration step. Additionally, first derivatives of the matrix with respect to parameters are needed at each step. Numerical accuracy of the blo ck-diagonalization depends on the separation sep(S, S ) of the diagonal blo cks in the Schur canonical form (calculated prior the blo ck-diagonalization) [15]. Instability o ccurs for very small values of sep(S, S ), which indicates that the spectra of S and S overlap under a very small perturbation of A0 . Thus, numerical instability signals that the chosen set of d eigenvalues should be changed such that the matrix S includes all "interacting" eigenvalues. The functions qi (p) are strongly nonlinear near the boundary of the surface d . The boundary corresponds to higher co dimension strata asso ciated with eigenvalues of higher multiplicity (or eigenvalues of the same multiplicity but with several Jordan blo cks). For example, the stratum 2 in Figure 1 is bounded by the singularities 3 and 4 . As a result, the convergence of Newton's metho d may be po or near the boundary of d . This instability signals that we should lo ok for eigenvalues with a more degenerate Jordan structure (e.g. higher multiplicity d). Analysis of the surface d very close to the boundary is still possible, but the higher precision arithmetics may be necessary. Figure 3 shows first iterations of Newton's pro cedure for different initial points in parameter space for matrix family (5.1) from Example 1. Solid arrows lo cate double eigenvalues (the stratum 2 ) and the dashed arrows correspond to triple eigenvalues (the stratum 3 ). One can see that the stratum 2 is well approximated when p0 is relatively far from the singularity (from the more degenerate stratum 3 ). For the left-most point p0 in Figure 3, the nearest point p 2 simply do es not exist (infimum of the distance p - p0 for p 2 corresponds to the origin p = 0 3 ). Note that, having the information on 12


Figure 3: One-step approximations of the nearest points of the strata 2 (solid arrows) and 3 (dashed arrows) near the cusp singularity. the stratum 3 , it is possible to determine lo cally the bifurcation diagram (describe the geometry of the cusp singularity in parameter space) [27, 35]. For the backward error analysis of numerical eigenvalue problems based on the study of the pseudo-spectrum we refer to [37]. We note that the classical numerical eigenvalue problem is ill-conditioned in the presence of multiple eigenvalues. The reason for that is the nonsmo othness of eigenvalues at multiple points giving rise to singular perturbation terms of order 1/d , where d is the size of Jordan blo ck [35]. On the contrary, in our problem we deal with the regular smo oth ob jects: the strata d and the versal deformation B(p).

7

Approximations based on diagonal decomp osition

In this section we consider approximations derived by using the diagonal decomposition of A0 . The diagonal decomposition is known to be ill-conditioned for nearly defective matrices. However, this way is easy to implement, while the very high accuracy may be not necessary. According to bifurcation theory for eigenvalues [35], the accuracy of the results based on the diagonal decomposition will be of order 1/d , where is the arithmetics precision. Another reason is theoretical. Bifurcation theory describes the collapse of a Jordan blo ck into simple eigenvalues [32, 35, 38]. Our approximations based on the diagonal decomposition solve the inverse problem: using simple (perturbed) eigenvalues and corresponding eigenvectors, we approximate the stratum d at which these eigenvalues coalesce. Let us assume that the matrix A0 is diagonalizable (its eigenvalues 1 , . . . , m are distinct). The right and left eigenvectors of A0 are determined by the equations A0 xi = i xi ,
yi A0 = i yi , yi xi = 1

(7.1)

with the last equality being the normalization condition. In Theorem 3.1 we take S = diag(1 , . . . , d ) = Y A0 X, X = [x1 , . . . , xd ], and Y = [y1 , . . . , yd ], where 1 , . . . , d are the eigenvalues coalescing at p d . In this 13


case expressions (3.5) take the form 1 q1 = pj d qi = pj
d d yi i=1

A xi , pj
i-1 yi

i=1

(i - q1 (p0 ))

A xi - trace(C pj

i-1

q1 ) - pj

i-1

trace(C
k =2

i-1

Ek 1 )

qk , pj

(7.2)

i = 2, . . . , d; j = 1, . . . , n. The interesting feature of these expressions is that they depend only on the simple eigenvalues 1 , . . . , d and their derivatives with respect to parameters at p0 [35]: i A = yi xi , pj pj i = 1, . . . , d . (7.3)

For example, for d = 2 we obtain the first-order approximation of the surface 2 in the form of one linear equation
n

+
j =1

2

y2

A A x2 - y1 x1 (pj - p0 ) = 0, j pj pj

=

2 - 1 . 2

(7.4)

Let us intro duce the gradient vectors i = ( i / p1 , . . . , i / pn ), i = 1, 2, with the derivatives i / pj given by expression (7.3). Then the solution of (7.4) approximating the vector p 2 nearest to p0 is found as p = p0 - 2 - 1 2 - 1
2

.

(7.5)

It is instructive to compare this result with the first-order approximation of the nearest p 2 , if we consider 1 and 2 as smo oth functions
n

i (p) = i +
j =1

i (pj - p0 ) + O ( p - p0 2 ), j pj

i = 1, 2 .

(7.6)

Using (7.6) in the equation 1 (p) = 2 (p) and neglecting higher order terms, we find
n

j =1

2 1 - pj pj

(p j - p 0 ) = 1 - 2 , j

(7.7)

which yields the nearest p as p = p0 - Comparing (7.8) with (7.5), functions, we find the correct dir in the distance to the stratum 2 of the bifurcation taking place at eigenvalues and eigenvectors [32, 2 - 1 2 - 1
2

2.

(7.8)

we see that considering simple eigenvalues as smo oth ection to the nearest point p 2 , but make a mistake overestimating it exactly twice. This is the consequence p 2 and resulting in O ( p - p0 1/2 ) perturbation of 35, 38]. 14


8

C o n c lu s io n

In the paper, we developed Newton's metho d for finding multiple eigenvalues with one Jordan blo ck in multiparameter matrix families. The metho d provides the nearest parameter vector with a matrix possessing an eigenvalue of given multiplicity. It also gives the generalized eigenvectors and describes the lo cal structure (tangent plane) of the stratum d . The motivation of the problem comes from applications, where matrices describe behavior of a system depending on several parameters. The whole matrix space has been considered as a particular case, when all entries of a matrix are independent parameters. Then the metho d provides an algorithm for solving the Wilkinson problem of finding the distance to the nearest degenerate matrix. Only multiple eigenvalues with one Jordan blo ck have been studied. Note that the versal deformation is not universal for multiple eigenvalues with several Jordan blo cks (the functions q1 (p), . . . , qd (p) are not uniquely determined by the matrix family) [3, 4]. This requires mo dification of the metho d. Analysis of this case is the topic for further investigation.

9
9.1

Ap p e n d i x
Pro of of Theorem 3.1
A0 U0 = U0 B0 , (9.1)

Taking equation (2.5) at p0 , we obtain

where U0 = U(p0 ) and B0 = B(p0 ). Comparing (9.1) with (3.1), we find that the matrix B0 is equivalent up to a change of basis to the matrix S. Then the equality (3.3) is obtained by equating the traces of the matrices B0 and S, where B0 has the form (2.5). Similarly, the equality (3.4) is obtained by equating the characteristic equations of the matrices B0 - q1 (p0 )I and S - q1 (p0 )I. The columns of the matrices X and U0 span the same invariant subspace of A0 . Hence, the matrices X and U0 are related by the expression U0 = XF, (9.2)

for some nonsingular d â d matrix F. Using (3.1) and (9.2) in (9.1), we find the relation SF = FB0 . (9.3) Taking derivative of equation (2.5) with respect to parameter pj at p0 , we obtain A
0

U U B A - B0 = U0 - U0 . pj pj pj pj

(9.4)

Let us multiply both sides of (9.4) by the matrix F-1 (S - q1 (p0 )I)i-1 Y and take the trace. Using expressions (3.1), (9.3), and the property trace(AB) = trace(BA), it is 15


straightforward to check that the left-hand side vanishes and we obtain the equation 0 = trace U0 B A - U0 F-1 (S - q1 (p0 )I) pj pj
i-1

Y



.

(9.5)

Substituting (9.2) into (9.5) and using equalities (3.1), (9.3), and B0 = q1 (p0 )I + C, we find B A trace Ci-1 - trace (S - q1 (p0 )I)i-1 Y X = 0. (9.6) pj pj Using (2.5), (3.6) and taking into account that trace(C trace(C
i-1 i-1

Ek1) = 0 for k > i, we obtain
i-1

)

q1 + pj

i

trace(C
k =2

i-1

Ek 1 )

qk = trace (S - q1 (p0 )I) pj

Y



A X. pj

(9.7)

Taking equation (9.7) for i = 1, . . . , d and using the equality trace(Ci-1 Ei1 ) = 1, we get the recurrent pro cedure (3.5) for calculation of derivatives of q1 (p), . . . , qd (p) at p0 .

9.2

Pro of of Theorem 3.2

At p0 d , we have q2 (p0 ) = · · · = qd (p0 ) = 0. From (2.5) it follows that B0 = J is the Jordan blo ck with the eigenvalue = q1 (p0 ) = trace S/d. By using (3.1), one can check that the vectors (3.10) satisfy the Jordan chain equations (2.1) for any vector k. Equations (3.11) provide the way of cho osing a particular value of the vector k. Since B0 = J , the versal deformation equation (2.5) becomes the Jordan chain equation (2.2) at p0 d . Hence, the columns of the matrix U0 are the generalized eigenvectors. Since the function q1 (p) and the matrix U(p) smo othly depend on parameters, the accuracy of the multiple eigenvalue and generalized eigenvectors has the same order as the accuracy of p0 d .

Acknowledgments
The author thanks A. P. Seyranian and Yu. M. Nechepurenko for fruitful discussions. This work was supported by the research grants RFBR 03-01-00161, CRDF-BRHE Y1-M-06-03, and the President of RF grant MK-3317.2004.1.

References
[1] E. Anderson, Z. Bai, C. Bischof, L. S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney and D. Sorensen. LAPACK Users' Guide (3rd edn). SIAM: Philadelphia, 1999. [2] I. E. Antoniou, M. Gadella, E. Hernandez, A. Jauregui, Yu. Melnikov, A. Mondragon and G. P. Pronko. Gamow vectors for barrier wells. Chaos, Solitons and Fractals 2001; 12:2719­2736. 16


[3] V. I. Arnold. Matrices depending on parameters. Russian Math. Surveys 1971; 26(2):29­43. [4] V. I. Arnold. Geometrical Methods in the Theory of Ordinary Differential Equations. Springer-Verlag: New York-Berlin, 1983. [5] M. V. Berry and M. R. Dennis. The optical singularities of birefringent dichroic chiral crystals. Proc. Roy. Soc. Lond. A 2003; 459:1261­1292. [6] J. V. Burke, A. S. Lewis and M. L. Overton. Optimal stability and eigenvalue multiplicity. Found. Comput. Math. 2001; 1(2):205­225. [7] J. W. Demmel. Computing stable eigendecompositions of matrices. Linear Algebra Appl. 1986; 79:163­193. [8] I. Dobson, J. Zhang, S. Greene, H. Engdahl and P. W. Sauer. Is strong mo dal resonance a precursor to power system oscillations? IEEE Transactions On Circuits And Systems I: Fundamental Theory And Applications 2001; 48:340­349. [9] A. Edelman, E. Elmroth and B. K°gstr¨m. A geometric approach to perturbation a o theory of matrices and matrix pencils. I. Versal deformations. SIAM J. Matrix Anal. Appl. 1997; 18(3):653­692. [10] A. Edelman and Y. Ma. Staircase failures explained by orthogonal versal forms. SIAM J. Matrix Anal. Appl. 2000; 21(3):1004­1025. [11] E. Elmroth, P. Johansson and B. K°gstr¨m. Bounds for the distance between nearby a o Jordan and Kronecker structures in closure hierarchy. In Numerical Methods and Algorithms XIV, Zapiski Nauchnykh Seminarov (Notes of Scientific Seminars of POMI) 2000; 268:24­48. [12] E. Elmroth, P. Johansson and B. K°gstr¨m. Computation and presentation of graphs a o displaying closure hierarchies of Jordan and Kronecker structures. Numerical Linear Algebra with Applications 2001; 8(6­7):381­399. [13] T. F. Fairgrieve. The application of singularity theory to the computation of Jordan canonical form. M.Sc. Thesis, Department of Computer Science, University of Toronto, 1986. [14] F. R. Gantmacher. The Theory of Matrices. AMS Chelsea Publishing: Providence, RI, 1998. [15] G. H. Golub and C. F. Van Loan. Matrix computations (3rd edn). Johns Hopkins University Press: Baltimore, 1996. [16] G. H. Golub and J. H. Wilkinson. Ill-conditioned eigensystems and the computation of the Jordan canonical form. SIAM Rev. 1976; 18(4):578­619.

17


[17] W. D. Heiss. Exceptional points ­ their universal o ccurrence and their physical significance. Czech. J. Phys. 2004; 54:1091­1099. [18] B. K°gstr¨m and A. Ruhe. An algorithm for numerical computation of the Jordan a o normal form of a complex matrix. ACM Trans. Math. Software 1980; 6(3):398­419. [19] B. K°gstr¨m and A. Ruhe. ALGORITHM 560: JNF, An algorithm for numerical a o computation of the Jordan normal form of a complex matrix [F2]. ACM Trans. Math. Software 1980; 6(3):437­443. [20] O. N. Kirillov. A theory of the destabilization paradox in non-conservative systems. Acta Mechanica 2005; to appear. [21] O. N. Kirillov and A. P. Seyranian. Optimization of stability of a flexible missile under follower thrust, in Proceedings of the 7th AIAA/USAF/NASA/ISSMO Symposium on Multidisciplinary Analysis and Optimization, St. Louis, MO, USA, 1998. AIAA Paper #98-4969. P. 2063­2073. [22] H. J. Korsch and S. Mossmann. Stark resonances for a double quantum well: crossing scenarios, exceptional points and geometric phases. J. Phys. A: Math. Gen. 2003; 36:2139­2153. [23] V. Kublanovskaya. On a metho d of solving the complete eigenvalue problem of a degenerate matrix. USSR Comput. Math. Math. Phys. 1966; 6(4):1­14. [24] O. Latinne, N. J. Kylstra, M. Dorr, J. Purvis, M. Terao-Dunseath, C. J. Jo chain, P. G. Burke and C. J. Noble. Laser-induced degeneracies involving autoionizing states in complex atoms. Phys. Rev. Lett. 1995; 74:46­49. [25] A. S. Lewis and M. L. Overton. Eigenvalue optimization. Acta Numerica 1996; 5:149­ 190. [26] R. A. Lippert and A. Edelman. The computation and sensitivity of double eigenvalues. In Advances in Computational Mathematics. Lecture Notes in Pure and Applied Mathematics. Vol. 202, Z. Chen, Y. Li, C. A. Micchelli, and Y. Xi, (eds). Dekker: New York, 1999; 353-393. [27] A. A. Mailybaev. Transformation of families of matrices to normal forms and its application to stability theory. SIAM J. Matrix Anal. Appl. 2000; 21(2):396­417. [28] A. A. Mailybaev. Transformation to versal deformations of matrices. Linear Algebra Appl. 2001; 337(1­3):87­108. [29] A. A. Mailybaev and A. P. Seiranyan. On the domains of stability of Hamiltonian systems. J. Appl. Math. Mech. 1999; 63(4):545­555. [30] A. A. Mailybaev and A. P. Seyranian. On singularities of a boundary of the stability domain. SIAM J. Matrix Anal. Appl. 2000; 21(1):106­128. 18


[31] MATLAB routines for computation of multiple eigenvalues and generalized eigenvectors for matrices dependent on parameters, on request by e-mail: mailybaev@imec.msu.ru [32] J. Moro, J. V. Burke and M. L. Overton. On the Lidskii-Vishik-Lyusternik perturbation theory for eigenvalues of matrices with arbitrary Jordan structure. SIAM J. Matrix Anal. Appl. 1997; 18(4):793­817. [33] A. Ruhe. An algorithm for numerical determination of the structure of a general matrix. BIT 1970; 10:196­216. [34] A. P. Seyranian, O. N. Kirillov and A. A. Mailybaev. Coupling of eigenvalues of complex matrices at diabolic and exceptional points. J. Phys. A: Math. Gen. 2005; 38, to appear. [35] A. P. Seyranian and A. A. Mailybaev. Multiparameter Stability Theory with Mechanical Applications. World Scientific: Singapore, 2003. [36] A. P. Seyranian and P. Pedersen. On interaction of eigenvalue branches in nonconservative multi-parameter problems. In Dynamics and Vibration of Time-Varying Systems and Structures: Conf. on Mech. Vibrat. and Noise. ASME: New York, 1993; 19­31. [37] E. Traviesas. Sur le d´ploiement du champ spectral d'une matrice. Ph.D. dissertation, e Universite Toulouse I, France, 2000. [38] M. I. Vishik and L. A. Lyusternik. The solution of some perturbation problems for matrices and selfadjoint or non-selfadjoint differential equations I. Russian Math. Surveys 1960; 15:1­74. [39] J. H. Wilkinson. The Algebraic Eigenvalue Problem. Clarendon Press: Oxford, 1965. [40] J. H. Wilkinson. Sensitivity of eigenvalues. Util. Math. 1984; 25:5­76. [41] J. H. Wilkinson. Sensitivity of eigenvalues II. Util. Math. 1986; 30:243­286.

19