This is a summary of my talk at the conference Celebrating the Centenary of James H. Wilkinson's Birth at the University of Manchester, May 29.

Jim Wilkinson knew that a matrix with badly conditioned simple eigenvalues must be close to a matrix with multiple eigenvalues. How can you find the nearest matrix where some of the eigenvalues have found their mates?

Much of this talk is my interpretation of the 2011 paper by Rafikul Alam, Shreemayee Bora, Ralph Byers and Michael Overton, "Characterization and construction of the nearest defective matrix via coalescence of pseudospectral components," in Linear Algebra and its Applications, vol. 435, pp. 494-513. Overton's accompanying code neardefmat is over 1000 lines of MATLAB® and is available from his web page.

Contents

A 3-by-3 Example

In last week's blog post, "An Eigenvalue Sensitivity Example", I resurrected this example from the 1980 User's Guide. For a moment, let's forget how it was constructed and start with

   A = [ -64 82 21; 144 -178 -46; -771 962 248]

A =
   -64    82    21
   144  -178   -46
  -771   962   248

The computed eigenvalues are

   format long
   lambda = eig(A)

lambda =
   3.000000000003868
   0.999999999998212
   1.999999999997978

The exact eigenvalues are 1, 2 and 3. We've lost about four figures. This is predicted by the eigenvalue condition numbers,

   format short
   kappa = condeig(A)

kappa =
  833.1092
  450.7228
  383.7564

These condition numbers come from the right eigenvector matrix

   [X,~] = eig(A)

X =
    0.0891    0.0735   -0.1089
   -0.1782   -0.1923    0.1634
    0.9800    0.9786   -0.9805

And the left eigenvector matrix

   Y = inv(X)

Y =
 -521.9612  628.5984  162.7621
  265.1820 -353.5760  -88.3940
 -257.0058  275.3634   73.4302

The rows of Y are left eigenvectors. They are normalized so that their dot product with the right eigenvectors is one.

   Y*X

ans =
    1.0000   -0.0000   -0.0000
   -0.0000    1.0000   -0.0000
   -0.0000    0.0000    1.0000

A function introduced in 2017b that computes the 2-norm of the columns of a matrix comes in handy at this point.

   kappa = (vecnorm(Y').*vecnorm(X))'

kappa =
  833.1092
  450.7228
  383.7564

EigTool is open MATLAB software for analyzing eigenvalues, pseudospectra, and related spectral properties of matrices. It was developed at Oxford from 1999 - 2002 by Thomas Wright under the direction of Nick Trefethen. It is now maintained by Mark Embree.

Here is the pseudospectrum of our example matrix.

eigtool(A)

You can see the three eigenvalues at 1, 2 and 3. The contours outline the regions where the eigenvalues can move when the matrix is perturbed by a specified amount.

This animated gif adjusts the region around the eigenvalues at 2 and 3 until the regions coalesce in a saddle point near 2.4.

To see more precisely where they coalesce, we use the function neardefmat, nearest defective matrix, written by Michael Overton, a professor at NYU. The outputs from neardefmat are B, the nearest defective matrix, and z, the resulting double eigenvalue. By default neardefmat prints some useful details about the computation.

   format long
   [B,z] = neardefmat(A)

neardefmat: lower bound on distance is 0.000299834 and upper bound is 0.00618718
neardefmat: finished, lowest saddle found is 2.41449 + 0 i with
 dist = 0.000375171, mu = 0, |u'*v| = 3.67335e-06
 singular vector residual norms are 1.8334e-13, 5.17298e-14
 this saddle point was found in search # 3
distance to nearest defective matrix B found is 0.000375171
neardefmat: two closest computed eigenvalues of B differ by 0.00269716
 and condition numbers of these eigenvalues are 272230, 272669
B =
   1.0e+02 *
  -0.639999774916363   0.819999582124016   0.210002343527353
   1.439999736641696  -1.779999511065701  -0.460002742035793
  -7.710000068576393   9.620000127314574   2.479999285995845
z =
   2.414492300748820

B is constructed from z, an SVD and a rank one perturbation of A.

   I = eye(3);
   [U,S,V] = svd(A - z*I);
   format short
   sigma = S(3,3)
   u = U(:,3)
   v = V(:,3)
   format long
   B = A - sigma*u*v'

sigma =
   3.7517e-04
u =
   -0.6373
    0.7457
    0.1942
v =
    0.0941
   -0.1748
    0.9801
B =
   1.0e+02 *
  -0.639999774916363   0.819999582124016   0.210002343527353
   1.439999736641696  -1.779999511065701  -0.460002742035793
  -7.710000068576393   9.620000127314574   2.479999285995845

This is essentially the same perturbation that I used to construct A in last week's blog post. The size of the perturbation is

   sigma = S(3,3)

   sigma = norm(B - A)

sigma =
     3.751705175569800e-04
sigma =
     3.751705175556183e-04

By construction, B has a double eigenvalue, but the eigenvalue's condition is infinite, so it is not computed accurately.

   eig(B)

ans =
   2.417189453582737
   2.414492295643361
   1.168318252152021

Frank(9)

Wilkinson was interested in the Frank family of matrices. Of course, we can find them in the gallery.

   help private/frank

 FRANK Frank matrix.
    F = GALLERY('FRANK',N, K) is the Frank matrix of order N.  It is
    upper Hessenberg with determinant 1. If K = 1, the elements are
    reflected about the anti-diagonal (1,N)--(N,1).  The eigenvalues
    of F may be obtained in terms of the zeros of the Hermite
    polynomials.  They are positive and occur in reciprocal pairs;
    thus if N is odd, 1 is an eigenvalue.  F has FLOOR(N/2)
    ill-conditioned eigenvalues---the smaller ones.

The 9x9 Frank matrix fits nicely on the page.

   F = gallery('frank',9)

F =
     9     8     7     6     5     4     3     2     1
     8     8     7     6     5     4     3     2     1
     0     7     7     6     5     4     3     2     1
     0     0     6     6     5     4     3     2     1
     0     0     0     5     5     4     3     2     1
     0     0     0     0     4     4     3     2     1
     0     0     0     0     0     3     3     2     1
     0     0     0     0     0     0     2     2     1
     0     0     0     0     0     0     0     1     1

The matrix is clearly far from symmetric. Since it is almost upper triangular you might think the eigenvalues are near the diagonal elements, but they are not. Here they are, sorted.

   format long
   [lambda,p] = sort(eig(F));
   lambda

lambda =
   0.044802721845067
   0.082015890201792
   0.162582729486328
   0.374121140912600
   0.999999999999966
   2.672931012564413
   6.150714797051293
  12.192759202150013
  22.320072505788492

And here are the corresponding condition numbers.

   format short e
   kappa = condeig(F);
   kappa = kappa(p)

kappa =
   1.4762e+04
   2.5134e+04
   1.1907e+04
   1.5996e+03
   6.8615e+01
   3.6210e+00
   1.5268e+00
   2.3903e+00
   1.9916e+00

The three smallest eigenvalues have the largest condition numbers. We expect the nearest defective matrix to bring two of them together.

Eigtool of the Frank matrix.

eigtool(F)

Focus on the three smallest eigenvalues. The red contour shows the regions containing the two that are coalescing.

Find the nearest defective matrix.

   format long
   [B,z] = neardefmat(F);

neardefmat: lower bound on distance is 3.4791e-07 and upper bound is 8.76938e-05
neardefmat: finished, lowest saddle found is 0.06104 + 0 i with
 dist = 5.01214e-07, mu = 0, |u'*v| = 6.09133e-10
 singular vector residual norms are 3.6435e-15, 6.12511e-16
 this saddle point was found in search # 1
distance to nearest defective matrix B found is 5.01214e-07
neardefmat: two closest computed eigenvalues of B differ by 1.7808e-06
 and condition numbers of these eigenvalues are 3.62946e+08, 3.62937e+08

The two smallest eigenvalues have merged, but they are hyper-sensitive.

   sort(eig(B))

ans =
   0.061039353359070
   0.061041134158948
   0.168054858422710
   0.373370423277969
   1.000016557229929
   2.672931171283966
   6.150714794349185
  12.192759202129578
  22.320072505788680

If we're not too strict on the tolerance, the eigenvector matrix X is rank deficient, implying the matrix B is defective.

   [X,~] = eig(B);
   rank(X,1.e-8)

ans =
     8

We can get a picture of the perturbation with the imagesc function. First, the matrix itself. You can see its upper Hessenberg (nearly upper triangular) structure.

imagesc(F)

Now, the picture of the perturbation shows that it is concentrated in the upper right-hand corner.

imagesc(F-B)

Wilkinson(9)

What perturbation of a symmetric matrix makes its eigenvalues coalesce? A great example was provided by Wilkinson himself.

   help wilkinson

   W = wilkinson(9)

 WILKINSON Wilkinson's eigenvalue test matrix.
    W = WILKINSON(n) is J. H. Wilkinson's eigenvalue test matrix, Wn+. It
    is a symmetric, tridiagonal matrix with pairs of nearly, but not
    exactly, equal eigenvalues. The most frequently used case is
    WILKINSON(21).
 
    W = WILKINSON(n,CLASSNAME) returns a matrix of class CLASSNAME, which
    can be either 'single' or 'double' (the default).
 
    Example:
 
    WILKINSON(7) is
 
           3  1  0  0  0  0  0
           1  2  1  0  0  0  0
           0  1  1  1  0  0  0
           0  0  1  0  1  0  0
           0  0  0  1  1  1  0
           0  0  0  0  1  2  1
           0  0  0  0  0  1  3

    Reference page in Doc Center
       doc wilkinson


W =
     4     1     0     0     0     0     0     0     0
     1     3     1     0     0     0     0     0     0
     0     1     2     1     0     0     0     0     0
     0     0     1     1     1     0     0     0     0
     0     0     0     1     0     1     0     0     0
     0     0     0     0     1     1     1     0     0
     0     0     0     0     0     1     2     1     0
     0     0     0     0     0     0     1     3     1
     0     0     0     0     0     0     0     1     4

W is symmetric and tridiagonal with ones on the off-diagonals. In order to have a multiple eigenvalue such a matrix has to have a pair of zeroes somewhere on the sub- and super-diagonal. But this matrix does not have any small off-diagonal elements. Nevertheless, the two largest eigenvalues are very close to each other.

   lambda = eig(W)

lambda =
  -1.125422415673319
   0.254718759825861
   0.952584219075215
   1.822717080887109
   2.178284739549981
   3.177282919112892
   3.247396472578982
   4.745281240174140
   4.747156984469142

Since the matrix is symmetric, the left and right eigenvectors are orthogonal and transposes of each other. All the eigenvalues are perfectly well conditioned. Where should we look for the coalescing eigenvalues?

   kappa = condeig(W)

kappa =
   1.000000000000000
   1.000000000000000
   1.000000000000000
   1.000000000000000
   1.000000000000000
   1.000000000000000
   1.000000000000000
   1.000000000000000
   1.000000000000000

Eigtool of Wilkinson(9). There are only seven circles for a 9x9 matrix -- the two circles on the right each contain two eigenvalues.

eigtool(W)

Let's zoom way in on the rightmost circle.

The -3.04 contours are just barely kissing. This is a different kind of saddle point -- a nonsmooth one. You will have to study the Alam et al paper to see why this is so important.

Let's see what neardefmat has to say about Wilkinson's matrix.

   [B,z] = neardefmat(W);

neardefmat: A is normal within tolerance
 distance to nearest defective matrix B is 0.000937872
neardefmat: two closest computed eigenvalues of B differ by 2.01677e-09
 and condition numbers of these eigenvalues are 465036, 465036

Sure enough, W is symmetric and hence normal. A rather different algorithm is required.

The perturbed eigenvalues have small imaginary parts.

   eig(B)

ans =
 -1.125422415673318 + 0.000000000000000i
  0.254718759825861 + 0.000000000000000i
  0.952584219075214 + 0.000000000000000i
  1.822717080887109 + 0.000000000000000i
  2.178284739549983 + 0.000000000000000i
  3.177282919112890 + 0.000000000000000i
  3.247396472578982 + 0.000000000000000i
  4.746219112321641 + 0.000000001008386i
  4.746219112321641 - 0.000000001008386i

Finally, here are images of the symmetric, tridiagonal matrix W and the very nonsymmetric perturbation W-B that brings two of its eigenvalues together.

Get the MATLAB code

Published with MATLAB® R2018b