# Pejorative Manifolds of Polynomials and Matrices, part 2

In an unpublished 1972 technical report "Conserving confluence curbs ill-condition," Velvel Kahan coined the descriptive term *pejorative manifold*. In case you don't use it in everyday conversation, *pejorative* means "expressing contempt or disapproval."

Velvel's report concerns polynomials with multiple roots, which are usually regarded with contempt because they are so badly conditioned. But Velvel's key observation is that although multiple roots are sensitive to arbitrary perturbation, they are insensitive to perturbations that preserve multiplicity.

Part 1 was about polynomials. This part is about matrix eigenvalues.

### Contents

#### The manifold

The pejorative manifold $\mathcal{M}$ is now the set of all 6-by-6 matrices which have an eigenvalue of multiplicity three at $\lambda$ = 3. These are severe restrictions, of course, and $\mathcal{M}$ is a tiny subset of the set of all matrices. But if we stay within $\mathcal{M}$, life is not nearly so harsh.

#### Two matrices

The Jordan Canonical Form of a matrix is bidiagonal, with eigenvalues on the diagonal and 1's and 0's on the super-diagonal. In our situation here, each eigenvalue with multiplicity m has a single m-by-m Jordan block with 1's on the superdiagonal. The Jordan blocks for distinct eigenvalues are separated by a zero on the super-diagonal.

Our first matrix has a 3-by-3 block with $\lambda$ = 3, then a 1-by-1 block with $\lambda$ = 2, and finally a 2-by-2 block with $\lambda$ = 1, so the diagonal is

d = [3 3 3 2 1 1];

and the superdiagonal is

j = [1 1 0 0 1];

Here it is

J1 = diag(j,1) A1 = diag(d) + J1

J1 = 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 A1 = 3 1 0 0 0 0 0 3 1 0 0 0 0 0 3 0 0 0 0 0 0 2 0 0 0 0 0 0 1 1 0 0 0 0 0 1

Our second matrix moves one super-diagonal element to swap the multiplicities of $\lambda$ = 2 and $\lambda$ = 1.

J2 = diag([1 1 0 1 0],1) A2 = diag([3 3 3 2 2 1]) + J2

J2 = 0 1 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 1 0 0 0 0 0 0 0 0 0 0 0 0 0 A2 = 3 1 0 0 0 0 0 3 1 0 0 0 0 0 3 0 0 0 0 0 0 2 1 0 0 0 0 0 2 0 0 0 0 0 0 1

These two matrices were constructed to have the two polynomials from my last post as characteristic polynomials. There is no need to compute the zeros, they are on the diagonals.

p1 = charpoly(A1,'s'); p2 = charpoly(A2,'s');

$$ p1 = s^6-13\,s^5+68\,s^4-182\,s^3+261\,s^2-189\,s+54 $$

$$ p2 = s^6-14\,s^5+80\,s^4-238\,s^3+387\,s^2-324\,s+108 $$

#### Convex combination

The convex linear combination puts the weights on the superdiagonal and the new eigenvalue of the diagonal.

```
format short
A = 1/3*A1 + 2/3*A2
```

A = 3.0000 1.0000 0 0 0 0 0 3.0000 1.0000 0 0 0 0 0 3.0000 0 0 0 0 0 0 2.0000 0.6667 0 0 0 0 0 1.6667 0.3333 0 0 0 0 0 1.0000

#### Characteristic polynomial

Let's check that the characteristic polynomial is the same as the third polynomial we had last time.

```
p3 = charpoly(A,'s');
```

$$ p3 = s^6-\frac{41\,s^5}{3}+76\,s^4-\frac{658\,s^3}{3}+345\,s^2-279\,s+90 $$

#### Plot

The plots of the three polynomials show how the triple root is more sensitive than either double root, which are, in turn, more sensitive than any of the simple roots.

plot_polys(p1,p2,p3)

#### Similarity transformation

A similarity transformation preserves, but disguises, the eigenvalues. Because it's handy, I'll use the HPL-AI matrix from my blog post last December HPL-AI Benchmark.

```
format short
M = HPLAI(6,-1)
```

M = 1.1667 0.1429 0.1250 0.1111 0.1000 0.0909 0.2000 1.1667 0.1429 0.1250 0.1111 0.1000 0.2500 0.2000 1.1667 0.1429 0.1250 0.1111 0.3333 0.2500 0.2000 1.1667 0.1429 0.1250 0.5000 0.3333 0.2500 0.2000 1.1667 0.1429 1.0000 0.5000 0.3333 0.2500 0.2000 1.1667

#### Yes ma'am

Here's how we do a similarity transformation in MATLAB.

B = M*A/M

B = 3.0610 1.1351 0.0899 -0.1695 -0.1178 -0.2053 0.0405 3.1519 1.1018 -0.1919 -0.1296 -0.2244 0.1360 0.2922 3.2144 -0.1402 -0.0745 -0.1867 0.1096 0.3808 0.2786 1.8527 0.6013 -0.1919 0.2839 0.6709 0.4447 -0.1222 1.6349 0.1467 1.5590 1.5424 0.7469 -0.1300 -0.0449 0.7517

#### Eigenvalues

What did this do to the eigenvalues?

```
format long
e = eig(B)
```

e = 1.000000000000000 + 0.000000000000000i 1.666666666666671 + 0.000000000000000i 1.999999999999999 + 0.000000000000000i 3.000006294572211 + 0.000000000000000i 2.999996852713897 + 0.000005451273553i 2.999996852713897 - 0.000005451273553i

The simple roots have hardly moved. The root with multiplicity three is perturbed by roughly the cube root of precision.

```
format short
e3 = e(4:6) - 3;
r3 = abs(e3)
circle(e3)
```

r3 = 1.0e-05 * 0.6295 0.6295 0.6295

#### Eigenvectors

What about the eigenvectors?

[V,~] = eig(B); imagV = imag(V) realV = real(V)

imagV = 1.0e-05 * 0 0 0 0 0 0 0 0 0 0 0.3705 -0.3705 0 0 0 0 0.0549 -0.0549 0 0 0 0 0.0678 -0.0678 0 0 0 0 0.0883 -0.0883 0 0 0 0 0.1225 -0.1225 realV = -0.0601 -0.0519 -0.0904 -0.6942 0.6942 0.6942 -0.0664 -0.0590 -0.1017 -0.1190 0.1190 0.1190 -0.0742 -0.0683 -0.1162 -0.1487 0.1487 0.1487 -0.3413 -0.9310 -0.9488 -0.1983 0.1983 0.1983 0.2883 0.3258 -0.1627 -0.2975 0.2975 0.2975 -0.8870 -0.1275 -0.2033 -0.5950 0.5950 0.5950

The last two vectors have small imaginary components and their real components are nearly the same as the fourth vector. So only the first four columns of `V` are good eigenvectors. We see `B`, and hence `A`, is *defective*. It does not have a full set of linearly independent eigenvectors.

#### Condition of eigenvalues

help condeig format short e kappa = condeig(B)

CONDEIG Condition number with respect to eigenvalues. CONDEIG(A) is a vector of condition numbers for the eigenvalues of A. These condition numbers are the reciprocals of the cosines of the angles between the left and right eigenvectors. [V,D,s] = CONDEIG(A) is equivalent to: [V,D] = EIG(A); s = CONDEIG(A); Large condition numbers imply that A is near a matrix with multiple eigenvalues. Class support for input A: float: double, single See also COND. Documentation for condeig doc condeig kappa = 1.6014e+00 2.9058e+00 2.9286e+00 1.3213e+10 1.3213e+10 1.3213e+10

Look carefully at those numbers -- the power of 10 on the first three is zero, while on the last three it is 10. This confirms that the first three eigenvalues are well conditioned, while the fourth is not.

#### Jordan Canonical Form

`A` is almost its own JCF, so this is no surprise.

A_jcf = jordan(A)

A_jcf = 3.0000e+00 1.0000e+00 0 0 0 0 0 3.0000e+00 1.0000e+00 0 0 0 0 0 3.0000e+00 0 0 0 0 0 0 1.0000e+00 0 0 0 0 0 0 1.6667e+00 0 0 0 0 0 0 2.0000e+00

But what about `B`? With exact computation, it would have the same JCF.

format short e B_jcf = jordan(B); real_B_jcf = real(B_jcf) imag_B_jcf = imag(B_jcf)

real_B_jcf = 1.0000e+00 0 0 0 0 0 0 1.6667e+00 0 0 0 0 0 0 2.0000e+00 0 0 0 0 0 0 3.0000e+00 0 0 0 0 0 0 3.0000e+00 0 0 0 0 0 0 3.0000e+00 imag_B_jcf = 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -2.7873e-06 0 0 0 0 0 0 2.7873e-06

The computed JCF is diagonal. This is, yet again, an example of the fact that the The Jordan Canonical Form Just Does't Compute. Also see Golub and Wilkinson.

#### And The Beat Goes On ...

I was just about to make that the end of the post when I tried this. I didn't realize it at the time, but the use of the HPL-AI matrix for the "random" similarity transformation has some welcome consequences. The elements of this matrix are the ratios of small integers.

M = sym(M)

M = [ 7/6, 1/7, 1/8, 1/9, 1/10, 1/11] [ 1/5, 7/6, 1/7, 1/8, 1/9, 1/10] [ 1/4, 1/5, 7/6, 1/7, 1/8, 1/9] [ 1/3, 1/4, 1/5, 7/6, 1/7, 1/8] [ 1/2, 1/3, 1/4, 1/5, 7/6, 1/7] [ 1, 1/2, 1/3, 1/4, 1/5, 7/6]

The elements of `A` are also the ratios of small integers.

sym(A)

ans = [ 3, 1, 0, 0, 0, 0] [ 0, 3, 1, 0, 0, 0] [ 0, 0, 3, 0, 0, 0] [ 0, 0, 0, 2, 2/3, 0] [ 0, 0, 0, 0, 5/3, 1/3] [ 0, 0, 0, 0, 0, 1]

The elements of `inv(M)` are ratios of large integers, but I don't have to show them because we're not going to invert `M`, even symbolically. Instead we use *forward slash* to compute the similarity transformation.

B = M*A/M;

The elements of `B` are also the ratios of large integers. Let's just look at the first column; the other columns are similar.

B1 = B(:,1)

B1 = 5261534240243927141/1718872313588352715 69679116377352174/1718872313588352715 46738873967260860/343774462717670543 2898457606578534/26444189439820811 97602262779214116/343774462717670543 2679724812276211392/1718872313588352715

#### JCF

With this exact symbolic computation, the characteristic polynomial of `B` is `p3`.

```
charpoly(B,'s')
```

ans = s^6 - (41*s^5)/3 + 76*s^4 - (658*s^3)/3 + 345*s^2 - 279*s + 90

As a result, the JCF of `B` is correct.

jordan(B)

ans = [ 1, 0, 0, 0, 0, 0] [ 0, 5/3, 0, 0, 0, 0] [ 0, 0, 2, 0, 0, 0] [ 0, 0, 0, 3, 1, 0] [ 0, 0, 0, 0, 3, 1] [ 0, 0, 0, 0, 0, 3]

#### Symbolic Eigenvectors

The third output argument of the symbolic version of `eig` is a vector whose length is the number of linearly independent eigenvectors. This is the sum of the *geometric multiplicities*. In this case, it is 4.

[V,E,k] = eig(B)

V = [ 11/27, 463/6831, 4/9, 7/6] [ 25/54, 31/414, 1/2, 1/5] [ 15/28, 485/5796, 4/7, 1/4] [ 460/63, 1115/2898, 14/3, 1/3] [ -23/9, -157/483, 4/5, 1/2] [ 1, 1, 1, 1] E = [ 5/3, 0, 0, 0, 0, 0] [ 0, 1, 0, 0, 0, 0] [ 0, 0, 2, 0, 0, 0] [ 0, 0, 0, 3, 0, 0] [ 0, 0, 0, 0, 3, 0] [ 0, 0, 0, 0, 0, 3] k = 1 2 3 4

#### Check

Verify that the eigenvectors work.

BV = B*V VE_k = V*E(k,k)

BV = [ 55/81, 463/6831, 8/9, 7/2] [ 125/162, 31/414, 1, 3/5] [ 25/28, 485/5796, 8/7, 3/4] [ 2300/189, 1115/2898, 28/3, 1] [ -115/27, -157/483, 8/5, 3/2] [ 5/3, 1, 2, 3] VE_k = [ 55/81, 463/6831, 8/9, 7/2] [ 125/162, 31/414, 1, 3/5] [ 25/28, 485/5796, 8/7, 3/4] [ 2300/189, 1115/2898, 28/3, 1] [ -115/27, -157/483, 8/5, 3/2] [ 5/3, 1, 2, 3]

#### TKP Preview

While preparing these two posts about "pejorative manifolds", I discovered some beautiful patterns of roundoff error generated by Triple Kronecker Products with very high multiplicity. That is going to be my next post. Here is a preview. tkp_preview.gif

close

## 댓글

댓글을 남기려면 링크 를 클릭하여 MathWorks 계정에 로그인하거나 계정을 새로 만드십시오.