Sensitivity and Computation of a Defective Eigenvalue

This website is set up to provide source codes and test results for the Algorithm PseudoEig developed in the paper by  Zhonggang Zeng

Matlab m-files for Algorithm PseudoEig:

PseudoEig.m,   peigmap.m,  peigjac.m, PEigRefine.m    (The NAClab package is required )

## Intuitive implementation of Algorithm PseudoEig             using NAClab

The NAClab package (Numerical Algebraic Computing laboratory) provides a quick and intuitive platform for implementing Algorithm PseudoEig.

For details of the algorithm and its theory, see the preprint  by  Zhonggang Zeng

A short sketch of the algorithm is a follows. NAClab provides an intuitive and straightforward platform to implement such an algorithm.  Just follow the steps below.

Step 2. Write a Matlab function for the mapping g and save it as peigmap.m , which is almost WYSIWYG.

Step 3.  Write a Matlab function for the Jacobian linear transformation and save it as peigjac.m, which is also almost WYSIWYG.

Step 4.  Write the main function PseudoEig and save it as PseudoEig.m,   where the Gauss-Newton iteration is calling the generic module GaussNewton.m in NAClab in a single line

[Z,res,lcond] = GaussNewton({@peigmap,{1,ones(n,k)},{A,C,S}}, @peigjac, {lambda0,X0},1);

By design, the generic module GaussNewton expects input items
--- the cell includes  (i)  the Matlab function handle for the mapping, i.e.  @peigmap
(ii)  the domain C x C^{n x k} represented by the cell, i.e. {1, ones(n,k)}
(iii)  the parameters for the mapping in the cell, i.e.  {A,C,S}
--- the Matlab function handle for the Jacobian linear transformation, i.e.   @peigjac
--- the initial iterate cell, i.e.  {lambda0, X0}
--- the optional tracking flag, 0 or 1,
The output Z is the cell whose first component is the pseudo eigenvalue and the second is the X.

Examples of using PseudoEig and PEigRefine
(Installation of NAClab is required to repeat the examples)
• Example 1:   Using PseudoEig, including finding the multiplicity support
• Example 2.   Refine the pseudo-eigenvalue

Example 1.  Using PseudoEig, including finding the multiplicity support

NAClab function "MatrixWithGivenJordanForm" is convenient in quick construction of testing matrices:

>> [A,U,V] = MatrixWithGivenJordanForm([2,3],{[4,3],[2 2]});
>> A

A =

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

Apply conventional eigen-finder for initial eigenvalues:

>> eig(A)

ans =

3.000000000000003 + 0.000000063443468i
3.000000000000003 - 0.000000063443468i
2.999999999999997 + 0.000000021680931i
2.999999999999997 - 0.000000021680931i
2.000131207566168 + 0.000131197500271i
2.000131207566168 - 0.000131197500271i
2.000003852182530 + 0.000006672920661i
2.000003852182530 - 0.000006672920661i
1.999992295634929 + 0.000000000000000i
1.999868792433838 + 0.000131217576475i
1.999868792433838 - 0.000131217576475i

Pick an approximate eigenvalue, say 1.999868792433838 + 0.000131217576475i, Find the geometric multiplicity by SVD:

>> svd(A-(1.999868792433838 + 0.000131217576475i)*eye(size(A)))

ans =

7.486664684066285
5.573841967969352
3.961701858109521
2.441112093740693
2.253848558674539
0.773876810987828
0.652336598832444
0.341709210710030
0.089726975408410
0.000000000002225
0.000000000000000

Clearly, the geometric multiplicity is 2.   We can then determine the Segre anchor k by trial and error.  Start with, say k=2:

>> [lambda,X,C,S,res,lcond] = PseudoEig(A,1.999868792433838 + 0.000131217576475i,2,2);
Step   1:     residual =   1.82e-08
Step   2:     residual =   6.38e-09
Step   3:     residual =   1.66e-09
>> res,lcond
res =
1.800832453591654e-11
lcond =
9.970824679378054e+08

Tiny residual and large condition number means k=2 underestimates the Segre anchor.  So increase it to k=3:

>> [lambda,X,C,S,res,lcond] = PseudoEig(A,1.999868792433838 + 0.000131217576475i,2,3);
Step   1:     residual =   2.02e-07
Step   2:     residual =   6.93e-14   ...

>> res, lcond
res =
2.775557561562891e-16
lcond =
4.162970862477419e+03

It appears that k=3 is the correct Segre anchor, the pseudo-eigenvalue

>> lambda
lambda =
2.000000000000000 - 0.000000000000000i

The pseodo-eigenvalue approximates the exact eigenvalue to the last digit.  To confirm the Segre anchor, try increasing it to k=4:

>> [lambda,X,C,S,res,lcond] = PseudoEig(A,1.999868792433838 + 0.000131217576475i,2,4);
Step   1:     residual =   4.35e-02
Step   2:     residual =   5.11e-05
Step   3:     residual =   7.20e-05
Step   4:     residual =   4.24e-05
Step   5:     residual =   4.14e-05
Step   6:     residual =   4.14e-05
Step   7:     residual =   4.13e-05
Step   8:     residual =   4.14e-05
Step   9:     residual =   4.14e-05
Step  10:     residual =   4.14e-05

The residual increased dramatically from 1e-16 to 1e-5, confirming k=4 is an overestimation of the Segre anchor.  Since the norm of the Moore-Penrose inverse of X is small, no refinement is needed:

>> norm(pinv(X))
ans =
2.649718399297023

Example 2.  Refine the pseudo-eigenvalue

>> A =   [ 2           1           0           0           0
0          -8           1           0           0
0           0           2           1           0
0           0           0           2           1
0      -10000        1000        -100          12];

>> [lambda,X,C,S,res,lcond] = PseudoEig(A,1.9998,1,5);
Step   1:     residual =   7.54e-07
Step   2:     residual =   7.62e-13
Step   3:     residual =   8.60e-14
Step   4:     residual =   1.20e-13
Step   5:     residual =   9.32e-14
Step   6:     residual =   8.96e-14
>> lambda

lambda =

1.999999999992827

The backward error is not small enough:

>> res*norm(pinv(X))
ans =
5.018090688429607e-08

Apply PEigRefine to refine the result:

>> lambda0 = lambda;  C0 = C; S0 = S; X0 = X;
>> [lambda,X,C,S,res,lcond] = PEigRefine(A,lambda0,X0,C0,S0)
Step   1:     residual =   1.82e-12
Step   2:     residual =   1.82e-12
Step   3:     residual =   1.90e-14
Step   4:     residual =   1.82e-12
Step   5:     residual =   1.82e-12
Step   6:     residual =   1.82e-12
>> lambda
lambda =
2

The refined backward error is now tiny:

>> res*norm(pinv(X))
ans =
1.901074022838812e-14

Special structure of the Jacobian matrix

The Matlab function  peigmat.m  rearranges the Jacobian matrix in to a near upper triangular form, making it possible to take advantage of the structure to improve efficiency.

Without loss of generality, assume the matrix is in Hessenberg form.  For a 100x100 matrix A, 100x3 matrix C,  3x3 S:

>> M = peigmat(lambda,X,A,C,S);
>> size(M)

ans =

159   151

>> spy(M) Example 4:  Contradicting sensitivity measures of the same eigenvalue

>> A = [1503/500,2,201/200,-1001/1000,-1/500,-1/1000,-1/1000,-1;
5,2,5,-1,-2,-1,-1,0;
-2503/500,-3,-601/200,2001/1000,1501/500,2001/1000,1/1000,2;
-6,-1,-6,3,5,3,0,1;
-5,-1,-5,1,6,3,0,1;
1,0,1,0,-1,1,0,0;
-4,-2,-4,1,3,2,2,2;
5,0,5,-1,-2,-1,-1,2]

>> e = eig(A)

e =
2.001415092944351 + 0.002122713121581i
2.001415092944351 - 0.002122713121581i
1.998908958350295 + 0.002124009818083i
1.998908958350295 - 0.002124009818083i
2.002684635776857 + 0.000000000000000i
1.997667261633844 + 0.000000000000000i
2.000000055770595 + 0.000000000000000i
1.999999944229407 + 0.000000000000000i

By PseudoEig

>> [lambda,X,C,S,res,lcond] = PseudoEig(A,e(1),2,2);
Step   1:     residual =   9.84e-06
Step   2:     residual =   9.82e-11
Step   3:     residual =   6.38e-16
Step   4:     residual =   5.96e-16
Step   5:     residual =   4.44e-16
Step   6:     residual =   5.15e-16
Step   7:     residual =   5.05e-16
Step   8:     residual =   9.33e-16
>> lambda, lcond

lambda =

2.000000000000000 - 0.000000000000000i

lcond =

3.304315987763452e+02

The 2x2 condition number is manageable at 3.3e+02, while the sensitivity measured by spectral projector norm is huge at  4.4e14

NAClab function NumericalJordanForm calculates the numerical Jordan Canonical Form

>> [J,X,e,s,t] = NumericalJordanForm(A,1e-9);
>> format short
>> J

J =

2.0028         0         0         0         0         0         0         0
0    2.0000   -0.8321         0         0         0         0         0
0         0    2.0000    0.2536         0         0         0         0
0         0         0    2.0000   -1.1073         0         0         0
0         0         0         0    2.0000    5.2368         0         0
0         0         0         0         0    2.0000         0         0
0         0         0         0         0         0    2.0000   -5.3852
0         0         0         0         0         0         0    2.0000

With a high overall condition number

>> cond(X)

ans =

5.0956e+13

even though the Jordan structure for the eigenvalue 2.00000000000 is not highly sensitive

>> t(2,2)

ans =

178.2683

If we use a loose error tolerance 1e-5, the numerical Jordan form becomes

>> [J,X,e,s,t] = NumericalJordanForm(A,1e-5);
>> single(J)

ans =

2.0001249   0.4962730           0           0           0           0           0           0
0   2.0001249   0.8362637           0           0           0           0           0
0           0   2.0001249  -0.3181596           0           0           0           0
0           0           0   2.0001249   1.1284051           0           0           0
0           0           0           0   2.0001249  18.9835510           0           0
0           0           0           0           0   2.0001249           0           0
0           0           0           0           0           0   2.0001249   5.3860307
0           0           0           0           0           0           0   2.0001249

The residual is within the error tolerance and condition number is not large:

>> t(1), t(2)
ans =
2.128714409525071e-06
ans =
1.590906220366132e+02

Example 1.  Identifying the multiplicity support:

>> A

A =

0   4   0  -4   0  -2   1   0   0  -1  -1  -1  -1   2   1   0   0  -1   0   0
0   3   3  -4   1   0   4   1  -1  -1   0   0   0   0   0  -2   0   0  -1  -1
1  -4   2  10  -3   1  -7  -2  -1   3   2   0   0  -1   0   3  -1   1   1   2
-1  -1   2   5  -2  -1  -5  -1  -1   2   0  -1   0   1   0   2  -1  -1  -1   1
-1  -2   2   1   1   1  -1   0  -2   1   0   0   0   1   0   0  -1  -1   0   0
1   4   1 -12   4   2  13   3   0  -4   0   0  -2  -1   1  -6   1   1   0  -3
-1  -1   1   5  -2   0  -4  -2   0   1  -1   0   0   1   0   4   0  -1  -1   2
1   2  -4   0   1   0   1   1   4  -2  -1   1   0  -1   0   1   2   1   0   1
0  -5   2  10  -5  -1 -10  -2  -1   6   3  -2   0   0   0   3  -3   0   1   2
1   1   1  -1   2   2   4   0   1  -1  -1   2   0  -1   0   0   2   1  -1   0
1  -1   0   2   1   2   1   0   1  -1   3   2   0  -1   0   0   1   1  -1   0
-1  -3   0   5  -1   2  -4  -1   0   1  -1   4   4   1  -2   2   0  -1   0   1
-2   0   0  -1   0   0   0   0   0   0  -1   0   3   2   0   0   0  -1   0   0
-3   4  -1  -4   0  -2   1   0   0  -1  -1  -1  -2   5   2   0   0  -1   1   0
-2   0   0  -1   0   0   0   0   0   0  -1   0   0   2   3   0   0  -1   0   0
0   0  -2   3   1   2  -1  -1   2  -2  -2   2   0   0   0   5   2   0   0   1
6   3  -6   3   6   4   7   0   7  -7   1   5  -2  -6   1   0   8   6  -1   0
0   2  -4  -4   1  -1   4   1   0  -1   0  -1  -1   0   1  -2   0   3   4  -1
1  -4  -1  11  -4   1  -8  -3  -1   3   2   0   0  -1   0   4  -1   1   4   3
0   0  -1   1  -2   0  -1  -2   1   0   0   0   0   0   0   1   0   0   0   4

Use the Matlab eig to find eigenvalue estimates:

>> diag(D)

ans =

3.001426431341224 + 0.000000000000000i
3.000624722484589 + 0.000453755506064i
3.000624722484589 - 0.000453755506064i
3.000440399461022 + 0.001356489469455i
3.000440399461022 - 0.001356489469455i
2.999761500834580 + 0.000734576671243i
2.999761500834580 - 0.000734576671243i
2.998846384868383 + 0.000837895621509i
2.998846384868383 - 0.000837895621509i
2.999227553361644 + 0.000000000000000i
1.999831389084806 + 0.000168679704578i
1.999831389084806 - 0.000168679704578i
2.000168610915185 + 0.000168542182678i
2.000168610915185 - 0.000168542182678i
1.999989850392488 + 0.000017578122305i
1.999989850392488 - 0.000017578122305i
1.999993977339806 + 0.000000000000000i
2.000020299214969 + 0.000000000000000i
2.000003011330126 + 0.000005216432308i
2.000003011330126 - 0.000005216432308i

The eigenvalue D(1,1) =  3.001426431341224 + 0.000000000000000i   appears to be defective since the condition number is high:

>> 1/abs(V(:,1)'*W(:,1))

ans =

3.096614239827069e+11

The geometric multiplicity is clearly two from the numerical nullity identified via singular values

>> s = svd(A-e(1)*eye(size(A)))

s =

40.024122556006361
22.819868881416550
11.340731455999819
9.397825324315559
6.260902607815194
4.358978843241264
3.295092277976436
2.208518179170260
2.014640549400185
1.237827273324948
0.661297941705164
0.546549612836940
0.316193846463476
0.194202964265260
0.108370106835163
0.070055683229182
0.054676439994994
0.036243966725033
0.000000000000002
0.000000000000001

To identify the Segre anchor k,  run PseudoEig with k = 1, 2, ...

>> for k = 1:6, [lambda,X,C,S,res,lcond] = PseudoEig(A,e(1),2,k); disp([k,res,lcond]); end

2.000000000000000e+00     1.467949463967694e-13     4.632744825087318e+14
3.000000000000000e+00     1.245977341450316e-10     1.515278877406559e+10
4.000000000000000e+00     2.365260962187022e-09     3.476672400527993e+09
5.000000000000000e+00     1.390184611834986e-15     5.791114059263334e+04
6.000000000000000e+00     2.622216751629970e-03     1.114112133754904e+05

The Segre anchor is  k = 5.

>> [lambda,X,C,S,res,lcond] = PseudoEig(A,e(1),2,5);
>> lambda

lambda =

3.000000000000086e+00

Example 2.  Refine the pseudo-eigenvalue by orthogonalization

By PseudoEig,  we get an approximate value of the eigenvalue

>> A

A =

2           1           0           0           0
0          -8           1           0           0
0           0           2           1           0
0           0           0           2           1
0      -10000        1000        -100          12

>> [lambda,X,C,S,res,lcond] = PseudoEig(A,1.999,1,5);
>> lambda

lambda =

1.999999999997139e+00

The backward error:

>> res*norm(pinv(X))

ans =

3.671544140185711e-08

Refine the pseudo-eigenvalue by orthogonalization:

>> lambda0 = lambda; X0 = X; C0 = C; S0 = S;
>> [lambda,X,C,S,res,lcond] = PEigRefine(A,lambda0,X0,C0,S0);
>> lambda

lambda =

2

>> res*norm(pinv(X))

ans =

2.874961099398239e-14

The backward error is reduced to 2.9e-14 from  3.7e-8

Example 3.  Refine pseudo-eigenvalue on a perturbed matrix

>> A

A =

2           1           0           0           0
0          -8           1           0           0
0           0           2           1           0
0           0           0           2           1
0      -10000        1000        -100          12

>> E

E =

-0.092000000000000  -0.653000000000000  -0.201000000000000  -0.416000000000000  -0.787000000000000
-0.135000000000000  -0.218000000000000   0.054000000000000  -0.136000000000000  -0.255000000000000
0.651000000000000   0.663000000000000  -0.166000000000000  -0.969000000000000  -0.603000000000000
-0.833000000000000   0.607000000000000   0.314000000000000   0.969000000000000  -0.020000000000000
-0.733000000000000  -0.879000000000000   0.256000000000000  -0.665000000000000  -0.321000000000000

>> B = A+1e-5*E

B =

1.0e+04 *

0.000199999908000   0.000099999347000  -0.000000000201000  -0.000000000416000  -0.000000000787000
-0.000000000135000  -0.000800000218000   0.000100000054000  -0.000000000136000  -0.000000000255000
0.000000000651000   0.000000000663000   0.000199999834000   0.000099999031000  -0.000000000603000
-0.000000000833000   0.000000000607000   0.000000000314000   0.000200000969000   0.000099999980000
-0.000000000733000  -1.000000000879000   0.100000000256000  -0.010000000665000   0.001199999679000

From the estimated eigenvalue 1.99

>> [U,T,V] = svd(B-1.99*eye(5));
N = V(:,end);
[lambda,X,C,S,res,lcond] = PseudoEig(B,1.99,1,5,N);
>> lambda

lambda =

2.004413315460551

Backward error

>> res*norm(pinv(X))

ans =

0.059004303654087

Refine:

>> lambda0 = lambda; X0 = X; C0 = C; S0 = S;
>> [lambda,X,C,S,res,lcond] = PEigRefine(B,lambda0,X0,C0,S0);
Step   1:     residual =   2.20e-02
Step   1:     residual =   1.89e-02
Step   2:     residual =   2.06e-06
Step   3:     residual =   2.06e-06
Step   4:     residual =   2.06e-06
Step   5:     residual =   2.06e-06
>> lambda

lambda =

1.999998778227778

With backward error:

>> res*norm(pinv(X))

ans =

9.603166323203087e-06