Jordan Canonical Form Zhonggang Zeng, Northeastern Illinois Univ. ( - - PowerPoint PPT Presentation

A Numerical Method

for Computing the

Jordan Canonical Form

Zhonggang Zeng, Northeastern Illinois Univ.

(Joint work with T. Y. Li, Michigan State Univ )

0 10 3 0 -1 -1 -4 0 0 -5 -5 0 1 0 0 -1 0 -5 -1 0 -3 -1 0 5 9 -1 3 -2 -1 1 1 -2 -2 1 -1 1 1 2 -1 -1 1 0 -1 1 3 1 7 2 -2 -11 1 0 6 -4 -3 6 0 5 -1 0 -3 -2 -1 0 0 0

• 1 1 5 2 3 1 -1 0 0 0 0 0 -1 0 1 2 0 0 1 0 -1 1
• 4 -2 -9 -2 6 19 -2 0 -8 8 6 -8 1 -7 1 -2 4 4 2 0 0 -1

0 -1 1 -1 1 2 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 9 -2 4 -3 3 3 1 -2 -2 1 0 1 2 1 -1 -1 1 0 -1 0 1 0 1 0 0 -2 0 3 4 0 0 3 0 2 0 0 -1 0 0 0 0 0 1 -4 -2 0 1 4 1 0 3 5 4 0 -2 0 0 1 0 3 1 0 1 1

• 1 1 -2 1 -1 3 -1 -1 -3 3 0 -3 0 -2 -1 0 1 0 0 0 0 0

5 2 6 2 -3 -16 1 0 12 -5 -1 12 0 9 -1 0 -5 -3 -2 0 0 0

• 1 4 0 1 -2 -4 -1 0 0 -5 -4 3 4 0 -1 -2 0 -3 -1 0 -1 -2

1 0 1 0 0 -2 0 0 2 0 0 2 3 2 0 0 -1 0 0 0 0 0 0 -1 4 -3 3 -1 1 1 0 0 0 0 -2 3 3 1 0 0 0 0 0 1 0 2 12 -1 2 -7 0 0 2 -4 -3 2 -3 2 4 6 -1 -2 0 0 -1 3

• 4 -1 -5 -2 2 12 -1 0 -7 4 3 -7 0 -6 1 3 4 2 1 0 0 0

0 11 8 1 -2 -12 -3 0 6 -9 -8 6 1 5 0 -1 0 -7 -2 0 -3 -1

• 2 0 7 -2 5 1 -1 1 -2 0 0 -2 0 -1 1 1 0 4 3 -1 -1 0

3 2 6 2 -2 -7 1 0 2 -5 -4 2 -2 2 0 3 -1 -3 1 2 0 2 5 -12 -10 2 -3 1 5 -1 0 6 6 0 0 0 -2 -1 0 6 0 3 5 0 4 -9 0 1 0 1 4 -1 0 4 4 0 -4 0 0 4 0 4 0 1 6 4 2 0 2 0 0 -3 0 0 3 0 0 3 0 3 0 0 -2 0 0 0 0 3 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3

JCF

RANMEP 2008, Jan 5, 2008, NCTS, Taiwan

Theorem: Every matrix corresponds to a Jordan Canonical Form

1 −

= X J X A

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ =

k k k k

J λ λ λ 1 1 O O

⎥ ⎥ ⎥ ⎥ ⎦ ⎤ ⎢ ⎢ ⎢ ⎢ ⎣ ⎡ =

l

J J J J O

2 1

1 However, it is extremely difficult to compute JCF with numerical computation

0 10 3 0 -1 -1 -4 0 0 -5 -5 0 1 0 0 -1 0 -5 -1 0 -3 -1 0 5 9 -1 3 -2 -1 1 1 -2 -2 1 -1 1 1 2 -1 -1 1 0 -1 1 3 1 7 2 -2 -11 1 0 6 -4 -3 6 0 5 -1 0 -3 -2 -1 0 0 0

• 1 1 5 2 3 1 -1 0 0 0 0 0 -1 0 1 2 0 0 1 0 -1 1
• 4 -2 -9 -2 6 19 -2 0 -8 8 6 -8 1 -7 1 -2 4 4 2 0 0 -1

0 -1 1 -1 1 2 1 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 1 9 -2 4 -3 3 3 1 -2 -2 1 0 1 2 1 -1 -1 1 0 -1 0 1 0 1 0 0 -2 0 3 4 0 0 3 0 2 0 0 -1 0 0 0 0 0 1 -4 -2 0 1 4 1 0 3 5 4 0 -2 0 0 1 0 3 1 0 1 1

• 1 1 -2 1 -1 3 -1 -1 -3 3 0 -3 0 -2 -1 0 1 0 0 0 0 0

5 2 6 2 -3 -16 1 0 12 -5 -1 12 0 9 -1 0 -5 -3 -2 0 0 0

• 1 4 0 1 -2 -4 -1 0 0 -5 -4 3 4 0 -1 -2 0 -3 -1 0 -1 -2

1 0 1 0 0 -2 0 0 2 0 0 2 3 2 0 0 -1 0 0 0 0 0 0 -1 4 -3 3 -1 1 1 0 0 0 0 -2 3 3 1 0 0 0 0 0 1 0 2 12 -1 2 -7 0 0 2 -4 -3 2 -3 2 4 6 -1 -2 0 0 -1 3

• 4 -1 -5 -2 2 12 -1 0 -7 4 3 -7 0 -6 1 3 4 2 1 0 0 0

0 11 8 1 -2 -12 -3 0 6 -9 -8 6 1 5 0 -1 0 -7 -2 0 -3 -1

• 2 0 7 -2 5 1 -1 1 -2 0 0 -2 0 -1 1 1 0 4 3 -1 -1 0

3 2 6 2 -2 -7 1 0 2 -5 -4 2 -2 2 0 3 -1 -3 1 2 0 2 5 -12 -10 2 -3 1 5 -1 0 6 6 0 0 0 -2 -1 0 6 0 3 5 0 4 -9 0 1 0 1 4 -1 0 4 4 0 -4 0 0 4 0 4 0 1 6 4 2 0 2 0 0 -3 0 0 3 0 0 3 0 3 0 0 -2 0 0 0 0 3

Problem : Jordan Canonical Form (JCF) 3

1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 3

JCF

) 10 (

15E

A

−

+ λ

{ }

3 ) ( = A λ

2

Computing JCF is an ill-posed problem

λ 1 λ λ 1 λ 1 λ λ 1 λ λ

A = X

E A A A ε + = → ~

Under arbitrary perturbation

λ1 λ λ2 λ5 λ3 λ7 λ6 λ8 λ4

Can we recover the lost eigenvalue and JCF

• f A using perturbed matrix à ?

3 λ 1 λ λ 1 λ 1 λ λ 1 λ λ

A X-1

λ8 λ6 λ7 λ5 λ4 λ3 λ2 λ1

X A ~ ~ =

1

~ − X

Background and motivation

• The new JCF algorithm is motivated by our recent works on solving

ill-posed problems in numerical polynomial algebra (approximate GCD, factorization, multiplicity structure, polynomial elimination, etc.)

• Numerical polynomial algebra and numerical algebraic geometry are fertile

application fields of numerical linear algebra

• JCF computing has direct application in numerical irreducible factorization
• f multivariate polynomials.

The preliminary implementation of the JCF algorithm is part of ApaLab (http://www.neiu.edu/~zzeng)

4

The pioneer work (Kublanovskaya-Ruhe):

• 2. Identify an eigenvalue cluster, and use the mean λ

as an approximation to a multiple eigenvalues

• 3. Use a sequence of SVD and unitary similarity transformations on A - λI

to obtain a local staircase form

• 1. Compute eigenvalues of A with Francis QR algorithm
• 4. Repeat on other multiple eigenvalues.

Reference: Kublanovskaya (1968) Ruhe (1970) Golub & Wilkinson (1976) Kågström & Ruhe (1980) Beelen & Van Dooren (1988) The main difficulty: accuracy of the multiple eigenvalue

5

>> J = find(abs(e-1)<0.2); >> sum(e(J))/20 ans = 1.00000000000000 >> K = find(abs(e-2)<0.2); >> sum(e(K))/15 ans = 2.00000000000000

If the cluster is properly identified, the center is usually be very good approximation to the multiple eigenvalue

6

The cluster mean can still be sensitive:

7

Can you solve (x-1.0 )100 = 0 x100-100 x99 +4950 x98 - 161700 x97+3921225x96 - ... - 100 x +1 = 0

8

A new development: Computing polynomial roots and multiplicities

“Computing multiple roots of inexact polynomials”,

• Z. Zeng, Math Comp, 2005

5 10 15 20

) 4 ( ) 3 ( ) 2 ( ) 1 ( − − − − x x x x

MultRoot results: The backward error: 6.16 x 10-16 Computed roots multiplicities 1.000000000000000 20 1.999999999999997 15 3.000000000000011 10 3.999999999999985 5

“Algorithm 835: MultRoot”,

• Z. Zeng, ACM TOMS, 2004

9

Question: Can we accurately compute multiple eigenvalues, and JCF?

(even if coefficients are approximate)

Why are ill-posed problems infinitely sensitive?

Plot of pejorative manifolds of degree 3 polynomials with multiple roots

• The solution structure is lost when the problem leaves the manifold due to an

arbitrary perturbation

• The problem may not be sensitive at all if the problem stays on the manifold,

unless it is near another pejorative manifold

10

• Problems with certain solution structure form a “pejorative manifold”
• W. Kahan’s observation (1972) on root-finding

{ }

) ( ) ( ) ( ) ( ) 1 , 1 , 1 (

1 1 1

γ β α γ β α ≠ ≠ − − − = = Π x x x x p

3

) 1 , 1 , 1 ( ) 2 , 1 ( ) 3 ( C = Π ⊂ Π ⊂ Π

Codimensions: 2 1 0 Stratification of pejorative manifolds

• f degree 3 polynomials

Pejorative manifold stratification of degree 4 polynomials:

) 4 ( Π ) 2 , 2 ( Π ) 3 , 1 ( Π ) 2 , 1 , 1 ( Π

4

) 1 , 1 , 1 , 1 ( C = Π

Codimensions: 3 2 1

{ }

β α β α ≠ − − = = Π

2 1

) ( ) ( ) ( ) 2 , 1 ( x x x p

{ }

C x x p ∈ − = = Π α α 3 ) ( ) ( ) 3 ( 11

If a polynomial p is near a pej. manifold, it is also near many pej. manifolds of lower codimesion

If p is a small perturbation from a pej. manifold ∏, then ∏ is the highest codimension manifold within a proper distance of p.

Pejorative manifolds (aka. bundles) of 4x4 matrices

e.g. {2,1} {1} is the structure of 2 eigenvalues in 3 Jordan blocks of sizes 2, 1 and 1 12

If a matrix A is near a bundle, it is also near many bundles of lower codimesion If A is a small perturbation from a bundle ∏, then ∏ is the highest codimension bundle within a proper distance of A.

{ }

) ( | matrices Rank * r A rank C A M r

n m n m r

= ∈ =

× ×

( )

) )( ( codim

• r

n r m M

n m r

− − =

×

{ }

)) , ( ( deg | ) , ( pairs al Polynomi *

,

r q p GCD q p P

n m r

= =

( )

r P

n m r

codim

• =

×

Geometry of ill-posed algebraic problems

n m n m n n m n

P P P

× × − ×

⊂ ⊂ ⊂

1

L

n m n n m n m

M M M

× × ×

⊂ ⊂ ⊂ L

1

Similar manifold stratification exists for problems such as polynomial factorization, multiple roots …

13

Objective of the approximate factorization: Given p(x), an approximate form of

k

m k m

z x z x x p ) ( ) ( ) ( ˆ

1

1

− − = L

find

k

m k m

z x z x x p ) ~ ( ) ~ ( ) ( ~

1

1

− − = L

) ( ˆ x p ) (x p ≈

Objective of the approximate JCF: Given a matrix A, an approximate form of

1

ˆ

−

= X J X A

find

1

~ ~ ~ ~

−

= X J X A A ≈

A A ~

14

1 codimsion = 2 n codimensio = 3 codimsion =

B

15

Illustration of pejorative manifolds

codimsion =

A

? ?

Problem A Problem B perturbation The “nearest” manifold may not be the answer The right manifold is of the highest codimension within a certain distance

Formulation of the approximate factorization:

Given p(x) and ε

> 0, its approximate factorization

k

m k m

z x z x x p ) ~ ( ) ~ ( ) ( ~

1

1

− − = L satisfies Backward nearness:

ε < − p p ~

Maximum codimension:

) ( ~

1 k

m m p L Π ∈

which is the of the highest codimension among all the factorization manifolds within ε

• f p

Minimum distance:

q p p p

k

m m q

− = −

Π ∈ ) (

1

min ~

L

The approximate factorization becomes a well-posed problem

p

p ˆ p ~

Approximate factorization of p = exact factorization of

p ~

approximating factorization of p

ˆ

ε

16

Continuity of the approximate solution:

Small perturbation will not change the maximum co-dimension manifold (if the given problem is near the manifold)

17

A “three-strikes” principle for formulating an approximate solution to an ill-posed algebraic problem:

• Backward nearness: The approximate solution is the exact solution of a nearby problem
• Maximum codimension: The approximate solution is the exact solution of a problem
• n the nearby pejorative manifold of the highest codimension.
• Minimum distance: The approximate solution is the exact solution of the nearest problem
• n the nearby pejorative manifold of the highest codimension.

Finding approximate solution is (likely) a well-posed problem Approximate solution is a generalization of exact solution.

18

after formulating the approximate solution of problem P within ε

P

The two-staged strategy

Stage II: Find/solve problem Q such that

R P Q P

R

− = −

Π ∈

min

Q Stage I: Find the pejorative manifold Π

• f the highest dimension s.t.

ε < Π) , (P dist Π

Exact solution of Q is the approximate solution of P within ε which approximates the solution of S where P is perturbed from S

19

by applying numerical linear algebra by the Gauss-Newton iteration

Univariate factorization: Stage I: Find the max-codimension pejorative manifold by applying univariate AGCD algorithm on (f, f’ )

n f x C f = > ∀ ∈ ∀ ) deg( , ], [ ε

1 m 1 m 1 1 m 1 m 1 m m 1

k 1 k 1 k 1

) ( ) ( ) ' , ( ) ( ) ( ) ( ) ( ' ) ( ) ( ) (

− − − −

− − ≈ ⇒ − − ≈ ⇒ − − ≈

k k k

z x z x f f AGCD x q z x z x x f z x z x x f L L L Q Stage II: solve the (overdetermined) polynomial system F(z1 ,…,zk )=f for a least squares solution (z1 ,…,zk ) by Gauss-Newton iteration (key theorem: The Jacobian is injective.)

) ( ) ( ) (

k 1

m m 1

• =

−

• −
• f

z z

k

L

(in the form of coefficient vectors)

20

Stage I: Identifying the multiplicity structure

p(x)= (x-1)5(x-2) 3(x-3) = [(x-1)4(x-2) 2]

[(x-1)(x-2)(x-3)]

p’(x)= [(x-1)4(x-2) 2]

[ (x-1)(x-2)+5(x -2)(x -3)+3(x -1)(x -3) ]

GCD(p,p’) = [(x-1)4(x-2) 2]

u0 (x) = [(x-1)4(x-2) 2]

[(x-1)(x-2)(x-3)]

u1 (x) = [(x-1)3(x-2) ]

[(x-1)(x-2)]

u2 (x) = [(x-1)2]

[(x-1)(x-2)]

u3 (x) = [(x-1)]

[(x-1)]

u4 (x) = [1]

[(x-1)]

distinct roots:

* * * * * * * * *

• multiplicities 5

3 1 u0 = p um =GCD(um-1 , um-1 ’)

21

Let ( x - z1 )

m1 (

x - z2 )

m2 ... (

x - zk )

mk =

xn + g1 (

z1 , ..., zk )

xn-1+...+gn-1 (

z1 , ..., zk )

x + gn (

z1 , ..., zk )

Then, p(x) = ( x - z1 )

m1 (

x - z2 )

m2 ... (

x - zk )

mk

<==> g1 (

z1 , ..., zk )

=a1 g2(

z1 , ..., zk )

=a2

... ... ...

gn (

z1 , ..., zk )

=an I.e. An over determined polynomial system

G(z) = a

(n > k)

(degree)

n k (number of distinct roots )

To minimize the distance from p(x)= xn + a1 xn-1+...+an-1 x + an to Π(m1 …mk )

22

Theorem: J(z) is of full rank <=> z1 ,…,zk are distinct. Theorem: J(z) is of full rank <=> z1 ,…,zk are distinct. ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎝ ⎛ = ) , , ( ) , , ( ) (

1 1 1 k n k

z z g z z g z G L M L

, The coefficient operator:

⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎜ ⎝ ⎛ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ =

k n n k

z g z g z g z g z J L M O M M O M L

1 1 1 1

) (

Its Jacobian:

Or the decomposition (

x - z1 )

m1 (

x - z2 )

m2 ... (

x - zk )

mk is irreducible

23

Theorem: The Gauss-Newton iteration locally converges

• at quadratic rate if the polynomial is exact
• at linear rate if the polynomial is inexact but close

The Gauss-Newton iteration z (i+1) =z(i)

• J(z

(i)

)+[ G(z (i) ) - a ],

i=0,1,2 ...

where J(.)+ is the pseudo-inverse of J(.)

24

5 10 15 20

) 4 ( ) 3 ( ) 2 ( ) 1 ( − − − − x x x x

For polynomial

with (inexact ) coefficients in machine precision

Stage I results: The backward error: 6.05 x 10-10 Computed roots multiplicities 1.000000000000353 20 2.000000000030904 15 3.000000000176196 10 4.000000000109542 5 Stage II results: The backward error: 6.16 x 10-16 Computed roots multiplicities 1.000000000000000 20 1.999999999999997 15 3.000000000000011 10 3.999999999999985 5

25

Formulation of the approximate Jordan Canonical Form:

Given matrix A ∈ Cn×n and ε

> 0, the approximate JCF

• f A = the exact JCF of

Backward nearness:

ε < − A A ~

Maximum codimension:

Π ∈ A ~

which is the of the highest codimension among All the matrix bundles within ε

• f A

Minimum distance:

B A A A

B

− = −

Π ∈

min ~

Conjecture: Computing the Jordan Canonical Form becomes a well-posed problem

A ~ A

A ˆ A ~

Approximate JCF of

A

= exact JCF of

A ~

approximating JCF of

A ˆ

ε

26

Segre and Weyr characteristics:

find X, J:

AX = XJ

Jordan form with Segre characteristics:

[ 3,2,2,1]

Staircase form with Weyr characteristics:

[ 4,3,1]

λ 1 λ λ 1 λ 1 λ λ 1 λ λ

J =

find U, S, λ

AU = U(λI+ S), UHU= I 3 2 2 1 4 3 1

Ferrer’s diagram

λI+ S=

λ λ λ λ λ λ λ λ + + + + + + + + + + + + + + + + + + + 27

A

The two-staged strategy

Stage II: Minimize the distance

B A A A

B

− = −

Π ∈

min ~

à Stage I: Find the higest codimension matrix bundle Π such that

ε < Π) , (A dist Π

Exact JCF of à is the approximate JCF of A within ε which approximates the JCF of  where A is perturbed from Â

28

by the Gauss-Newton iteration

find U, S, λ

AU = U(λI+S) UHU=I λI+ S=

λ λ λ λ λ λ λ λ + + + + + + + + + + + + + + + + + + +

Stage II: assume the Weyr characteristics is known, or computed at Stage I: for ( λ, U, S ) AU- U(λI+S) = O [c1 , …, cj ]Huj – ej = 0, j=1,…,m

leads to overdetermined quadratic system (with ci , bj at random)

(1)

with a staircase structure constraint

Φ ∈ = ) , ( , j i u b

i H j

29

that defines the bundle having an eigenvalue corresponds to the given Weyr characteristics.

Its least squares solution minimizes the distance to the bundle

Theorem: The Jacobian J(λ,U,S) of F(λ,U,S) is injective

• Eigenvalues are no longer flying!
• Gauss-Newton iteration converges locally

AU- U(λI+S) = O [c1 , …, cj ]Huj – ej = 0, j=1,…,m

Φ ∈ = ) , ( , j i u b

i H j

) , , ( = S U F λ

• Condition number:

∞ < <

+

) , , ( S U J λ

The Stage II problem is well-posed, even well-conditioned

30

Gauss-Newton iteration on system

) , , ( = S Y F λ

Example: A 50x50 matrix with eigenvalues 1, 2, 3 with multiplicity 20, 15, 5 respectively

31

SGMIN: Lippert & Edelman 2000 EigTrip: Zeng & Li, Stage II code

32

eigenvalue Segre characteristics λ1 5 3 2 2 λ2 3 3 1 λ3 2

• 10 6 3 2

Minimal polynomials:

2 1 4 1 2 2 1 3 3 2 3 1 2 2 3 3 2 5 1 1

) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ − = − − = − − = − − − = p p p p

Let v be a coefficient vector of p1 . Then for a random x

[ ]

, , ,

10

= v x A Ax x L

A rank/kernel problem

33

Stage I: Compute the Weyr/Segre characteristics

QH AQ =

A

=

1

q

…

k

q

1

q

…

1 + k

q

H

{ }

{ }

k k

q q span q A Aq q span , , , , ,

1 1 1 1 1

L L =

−

{ }

[ ] ( ) ( )

QH q Ran q q A q Ran q A Aq q span

k k

, , , , , , ,

1 1 1 1 1 1

= = L L

= H

1

Q Ran

Consider the Hessenberg reduction

[ ]

, , , = v x A Ax x

k

L

H is rank-deficient

1 34

We can further refine the Hessenberg reduction by Gauss-Newton iteration on

⎪ ⎩ ⎪ ⎨ ⎧ = − = − = −

1

x q O I Q Q O QH AQ

H

A

1

q

…

k

q

1

q

…

k

q

H

=

A Q

… … … …

H1 H2 Hm

QH

35

By a Hessenberg reduction,

with A1 = A, and a random x being the first column of Q,

Rank/kernel leads to p1 (λ)

(the 1st minimal polynomial)

leads to p2 (λ), p3 (λ), … recursively

36

Minimal polynomials

2 1 4 1 2 2 1 3 3 2 3 1 2 2 3 3 2 5 1 1

) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( ) ( λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ λ − = − − = − − = − − − = p p p p

approximate factorization with approximate data Rank-revealing GCD Root-finding JCF structure

37

Example: 100x100 matrix A with multiple eigenvalues 1, -1, 2, -2 50 simple eigenvalues: random

38

A(t) =

JNF: Kågström & Ruhe (1980) NumJCF: Zeng & Li

39

1 −

⎥ ⎦ ⎤ ⎢ ⎣ ⎡ = X B J X A

20 80

) 2 ( ) 2 ( ) 2 ( ) 1 ( ) 1 ( ) 1 ( ) 1 (

2 2 4 1 3 4 5

J J J J J J J J ⊕ ⊕ ⊕ ⊕ ⊕ ⊕ =

X and B are 100x100 and 80x80 random matrices, respectively

On 1000 such matrices, we run JNF and NumJCF twice each

40

Concluding remarks

• Computing the approx. JCF as formulated above appears to be a

well-posed problem, and numerical computation is possible

• A new algorithm is proposed, along with a condition number.
• Limit exists for computing JCF with perturbed data
• lots of work ahead