next up previous
Next: Remarks on the Up: Eigenvalue Problems. Previous: Symbolic Computation With

An Example Due to Kenton Yee.

Near the end of 1990, Kenton Yee of Brookhaven National Labs posted a question about the following matrix to the USENET newsgroup sci.math.symbolic. He wanted to know the eigenvectors of this matrix, so he could guess how to proceed with a similar but larger (32 by 32) matrix.

 

Naive use of Maple or any other CAS simply leads to difficulties. For example, Gerald Edgar pointed out in that series of postings that Maple 4.2 could correctly compute the eigenvalues in terms of radicals. They are 1, with multiplicity 3, and -1 again with multiplicity 3, and

Further, Maple could correctly compute , but it couldn't correctly calculate the nullspace in version 4.2. Maple V at least gets a nontrivial nullspace, but the answer is too complicated to understand:

> nullspace(B1);

                      1/2        2   1/2         5        6         3   1/2
{[ 2 (4 + 104 e + 3 %1    + 498 e  %1    + 2408 e  + 456 e  + 1351 e  %1

               4   1/2     3/2          1/2        2         3         4
        + 896 e  %1    + %1    + 67 e %1    + 904 e  + 2984 e  + 3484 e

               5   1/2        6   1/2     3/2          3/2  2       7   1/2
        + 693 e  %1    + 210 e  %1    + %1    e - 44 %1    e  + 25 e  %1

              3   3/2    8   1/2    4   3/2       7    /         2
        - 13 e  %1    + e  %1    - e  %1    + 24 e )  /  (e %4 %3  %2),
                                                     /

              1/2        2   1/2     5/2          5          6         3   1/2
1/2 (16 e - %1    - 672 e  %1    + %1    + 13120 e  + 10448 e  - 2754 e  %1

             4   1/2          1/2        2         3          4         5   1/2
     - 1938 e  %1    - 50 e %1    + 400 e  + 3312 e  + 10448 e  - 3054 e  %1

             6   1/2        3/2           3/2  2        7   1/2       3   3/2
     - 1048 e  %1    + 42 %1    e + 308 %1    e  - 126 e  %1    + 70 e  %1

          8   1/2      4   3/2         7        8       9    /         2
     - 5 e  %1    + 4 e  %1    + 3312 e  + 400 e  + 16 e )  /  (e %4 %3  %2),
                                                           /

            1/2       2   1/2        2        3            4       5    6
4 (20 e + %1    + 42 e  %1    + 113 e  + 164 e  + 1 + 113 e  + 20 e  + e

           3   1/2    4   1/2          1/2
     + 14 e  %1    + e  %1    + 14 e %1   )/(%4 %3 %2),                  1,

             1/2       2   1/2       2        3       4      5       3   1/2
1/2 (4 e + %1    + 38 e  %1    + 52 e  + 104 e  + 52 e  + 4 e  + 16 e  %1

        4   1/2     3/2          1/2
     + e  %1    - %1    + 16 e %1   )/(e %4 %3),

            1/2       2   1/2        2        3            4       5    6
4 (20 e + %1    + 42 e  %1    + 113 e  + 164 e  + 1 + 113 e  + 20 e  + e

           3   1/2    4   1/2          1/2
     + 14 e  %1    + e  %1    + 14 e %1   )/(%4 %3 %2),

            1/2        2   1/2        2        3            4       5      6
(68 e + 5 %1    + 106 e  %1    + 308 e  + 392 e  + 4 + 308 e  + 68 e  + 4 e

           3   1/2      4   1/2     3/2          1/2    /       2
     + 56 e  %1    + 5 e  %1    - %1    + 56 e %1   )  /  (%2 %3 ),
                                                      /

   2                 2       1/2    3       1/2
  e  (29 e + 5 + 13 e  + 7 %1    + e  + e %1   )
2 ---------------------------------------------- ]}
                       %3 %2

                            2              3        4
%1 :=                   10 e  + 12 e + 12 e  + 1 + e

                                       2     1/2
%2 :=                       8 e + 1 + e  + %1

                                       2     1/2
%3 :=                       4 e + 1 + e  + %1

                                        2     1/2
%4 :=                       10 e + 1 + e  + %1

Dealing with nested radicals is difficult. So it is not clear that the above nullspace is correct for all (or any) values of e. The following session uses the Row Echelon Decomposition and Maple's RootOf construct to explore this problem completely.

The following sequence of commands illustrates the use of the Row Echelon Decomposition [36] to solve the posted problem. The new feature of this solution is that special cases, which depend on the value of the parameter e, are automatically detected.

> with(linalg):
Warning: new definition for   trace
> with(share):
> readlib(Echelon):
Defines the routines RowEchelon, BackSubstitute, and some auxiliary routines.
> read `mat.mpl`;
> mat := array(1..8,1..8,
>   [[e^(-1), e^(-1), e^(-1), e^(-1), e^(-1), e^(-1), e^(-1), 0],
>    [1, 1, 1, 1, 1, 1, 0, 1], [1, 1, 1, 1, 1, 0, 1, 1],
>    [1, 1, 1, 1, 0, 1, 1, 1], [1, 1, 1, 0, 1, 1, 1, 1],
>    [1, 1, 0, 1, 1, 1, 1, 1], [1, 0, 1, 1, 1, 1, 1, 1],
>    [0, e, e, e, e, e, e, e]]);
                       [ 1/e  1/e  1/e  1/e  1/e  1/e  1/e  0 ]
                       [                                      ]
                       [  1    1    1    1    1    1    0   1 ]
                       [                                      ]
                       [  1    1    1    1    1    0    1   1 ]
                       [                                      ]
                       [  1    1    1    1    0    1    1   1 ]
                mat := [                                      ]
                       [  1    1    1    0    1    1    1   1 ]
                       [                                      ]
                       [  1    1    0    1    1    1    1   1 ]
                       [                                      ]
                       [  1    0    1    1    1    1    1   1 ]
                       [                                      ]
                       [  0    e    e    e    e    e    e   e ]


> A := evalm(mat - lambda*Idn(8));

              A := [1/e - lambda, 1/e, 1/e, 1/e, 1/e, 1/e, 1/e, 0]

                          [1, 1 - lambda, 1, 1, 1, 1, 0, 1]

                          [1, 1, 1 - lambda, 1, 1, 0, 1, 1]

                          [1, 1, 1, 1 - lambda, 0, 1, 1, 1]

                          [1, 1, 1, 0, 1 - lambda, 1, 1, 1]

                          [1, 1, 0, 1, 1, 1 - lambda, 1, 1]

                          [1, 0, 1, 1, 1, 1, 1 - lambda, 1]

                          [0, e, e, e, e, e, e, e - lambda]

> RowEchelon(A,'D','P','L','U'):
We find the eigenvalues by factoring the determinant of .
> chrp := factor(D);

chrp :=

            3        2             2                                          3
(lambda + 1)  (lambda  e - lambda e  - 6 lambda e - lambda + 7 e) (lambda - 1)
-------------------------------------------------------------------------------
                                       e
If e is zero, the matrix itself is not defined. This special case is not pursued further.

We see that 1 is a root of multiplicity 3, as is -1. Also, we have a quadratic factor which depends on e, so only two of the eigenvalues depend on e.

We find the eigenvectors for the eigenvalue , first.

> R1 := RowEchelon(subs(lambda=-1,op(A)),'D1');

                         [ 1  0  0  0  0  0   0     - 1/e   ]
                         [                                  ]
                         [ 0  1  0  0  0  0   -1      0     ]
                         [                                  ]
                         [ 0  0  1  0  0  -1  0       0     ]
                         [                                  ]
                         [                            e + 1 ]
                         [ 0  0  0  1  0  1   1   1/2 ----- ]
                         [                              e   ]
                   R1 := [                                  ]
                         [                            e + 1 ]
                         [ 0  0  0  0  1  1   1   1/2 ----- ]
                         [                              e   ]
                         [                                  ]
                         [ 0  0  0  0  0  0   0       0     ]
                         [                                  ]
                         [ 0  0  0  0  0  0   0       0     ]
                         [                                  ]
                         [ 0  0  0  0  0  0   0       0     ]
Notice that the rank of R is 5, so there are 3 eigenvectors associated with the eigenvalue -1.

We find the eigenvectors (all at once) by using a routine BackSubstitute, which allows for free parameters, which we will call T.i for i = 1,2,3,... (for no particular reason).

> ev1 := BackSubstitute(R1,'T');
                             T1              2 T3 e + 2 T2 e + T1 e + T1
       ev1 := array(1 .. 8,[----,T2,T3,- 1/2 ---------------------------,
                             e                            e

                 2 T3 e + 2 T2 e + T1 e + T1
           - 1/2 ---------------------------,T3,T2,T1])
                              e
Now check for special cases -- we look at the determinant D computed by RowEchelon. This is supposed to be nonzero. If it is zero for some value of the parameters, we have to re-do the calculation.
> normal(D1);
                                       4
Nonzero, independent of e. So these eigenvectors are ok, aside from the obvious problem at e=0. Now we look at the eigenvectors for lambda = +1.
> R2 := RowEchelon(subs(lambda=1,op(A)),'D2');

                              [ 1  0  0  0  0  0  0  0 ]
                              [                        ]
                              [ 0  1  0  0  0  0  1  0 ]
                              [                        ]
                              [ 0  0  1  0  0  1  0  0 ]
                              [                        ]
                              [ 0  0  0  1  1  0  0  0 ]
                        R2 := [                        ]
                              [ 0  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  0  0  0  0 ]

> ev2 := BackSubstitute(R2,'F');
               ev2 := array(1 .. 8,[0,- F2,- F3,- F4,F4,F3,F2,0])
Now again, check:
> normal(D2);
                                   2
                                  e  - 2 e + 1
                                  ------------
                                       e
Zero if and only if e = 1. This means we must repeat our Row Echelon Decomposition calculation, in the special case e=1.
> R2a := RowEchelon(subs([lambda=1,e=1],op(A)),'D2a');

                              [ 1  0  0  0  0  0  0  1 ]
                              [                        ]
                              [ 0  1  0  0  0  0  1  0 ]
                              [                        ]
                              [ 0  0  1  0  0  1  0  0 ]
                              [                        ]
                              [ 0  0  0  1  1  0  0  0 ]
                       R2a := [                        ]
                              [ 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  0  0  0  0 ]

> ev2a := BackSubstitute(R2a,'Q');
            ev2a := array(1 .. 8,[- Q1,- Q2,- Q3,- Q4,Q4,Q3,Q2,Q1])
Rank R2a = 4, so the eigenvector will have multiplicity 4. There are no free parameters left, so we know D2a <> 0. Now the hard case.
> quad_factor := normal(e*chrp/(lambda-1)^3/(lambda+1)^3);
                             2             2
        quad_factor := lambda  e - lambda e  - 6 lambda e - lambda + 7 e

> alias(alpha = RootOf(quad_factor,lambda));
Rather than using square roots and doing two calculations, we use SPM_quotRootOf" and do just one eigenvector calculation.

Note that SPM_quotalpha" here stands for either eigenvalue. It refers to one of them, not both. To put this rather slippery idea more clearly, think of it as telling Maple ``Ok, let be one of the roots of that quadratic. I know which one, but I am not telling you (yet). As far as you are concerned, the only simplifications you can do with that root are those valid for both roots."

This is an example of a `generalized representation of a function'. Given e, and a program for finding square roots efficiently and stably, this is a perfectly good description of each root.

Notice that working with this is much more difficult than the previous calculations. In general, recognizing when a functional expression is zero is an intractable problem. Working with algebraic numbers is not much easier. Here is an algebraic function of e, and Maple's (and other CAS') normalization facilities (= zero-recognition facilities) are often not strong enough to work on such problems. Certainly, the default normalizer used in RowEchelon (= Maple's SPM_quotnormal" function) is not strong enough.

This problem has an `echo' in the numerical analysis world. I think this is an example of the principle of `conservation of difficulty'. Given numerical values for e, then the decision of what the rank of the matrix is is theoretically easier, but in the presence of noise (roundoff error or uncertainty in e) it can be just as difficult. Matlab uses the SVD to decide numerical rank, and this procedure is robust even in the presence of roundoff or uncertainty.

Here, a sophisticated user of Maple would expect that the appropriate normalizer to use is SPM_quotevala@Normal". However, the algebraic function case was not yet implemented in Maple V, where this example was created. So we make do with a weaker normalizer, the combination of Maple's command SPM_quotsimplify", which knows about RootOf and can simplify things like and to expressions linear in , and SPM_quotnormal", which can handle rational functions.

> Normalizer := normal@simplify:

> R3 := RowEchelon(subs(lambda=alpha,op(A)),'D3');

               [                           2                         ]
               [                          e  - alpha e + 7 e - alpha ]
               [ 1  0  0  0  0  0  0  1/6 -------------------------- ]
               [                                       2             ]
               [                                      e              ]
               [                                                     ]
               [                                  e - alpha          ]
               [ 0  1  0  0  0  0  0          1/6 ---------          ]
               [                                      e              ]
               [                                                     ]
               [                                  e - alpha          ]
               [ 0  0  1  0  0  0  0          1/6 ---------          ]
               [                                      e              ]
               [                                                     ]
               [                                  e - alpha          ]
         R3 := [ 0  0  0  1  0  0  0          1/6 ---------          ]
               [                                      e              ]
               [                                                     ]
               [                                  e - alpha          ]
               [ 0  0  0  0  1  0  0          1/6 ---------          ]
               [                                      e              ]
               [                                                     ]
               [                                  e - alpha          ]
               [ 0  0  0  0  0  1  0          1/6 ---------          ]
               [                                      e              ]
               [                                                     ]
               [                                  e - alpha          ]
               [ 0  0  0  0  0  0  1          1/6 ---------          ]
               [                                      e              ]
               [                                                     ]
               [ 0  0  0  0  0  0  0                0                ]
We have succeeded (in a few minutes CPU time) in showing that the rank of is 7, for either . This handles both remaining cases at once. Note that the eigenvectors, found below, depend on e and depend, naturally, on the eigenvalue.
> ev3 := BackSubstitute(R3,'Q');
                             2
                           (e  - alpha e + 7 e - alpha) Q1       (e - alpha) Q1
ev3 := array(1 .. 8,[- 1/6 -------------------------------,- 1/6 --------------
                                          2                            e
                                         e

           (e - alpha) Q1       (e - alpha) Q1       (e - alpha) Q1
    ,- 1/6 --------------,- 1/6 --------------,- 1/6 --------------,
                 e                    e                    e

          (e - alpha) Q1       (e - alpha) Q1
    - 1/6 --------------,- 1/6 --------------,Q1])
                e                    e
Now we have to check that D3 is not zero.
> normal(simplify(D3));

      8            7       7              6        6               5
  - (e  alpha - 7 e  + 30 e  alpha - 168 e  + 333 e  alpha + 1692 e  alpha

               5         4         4               3               3         2
       - 1365 e  - 4368 e  + 4023 e  alpha - 5145 e  + 4470 alpha e  - 2520 e

                     2                                     /  2
       + 2663 alpha e  - 251 e + 576 alpha e + 36 alpha)  /  e
                                                         /
Now when is this zero? We can try SPM_quotsolve" directly, but it would be nice if we could find a polynomial expression equivalent to this, because working with nested roots is horrible.
> readlib(isolate)(",alpha);
  alpha =

              7        6         5         4         3         2
           7 e  + 168 e  + 1365 e  + 4368 e  + 5145 e  + 2520 e  + 251 e
      ------------------------------------------------------------------------
       8       7        6         5         4         3         2
      e  + 30 e  + 333 e  + 1692 e  + 4023 e  + 4470 e  + 2663 e  + 576 e + 36
is a function of e. If happens to be equal to the right-hand side above, then D3 is zero. Also, if D3 = 0, then happens to equal the right-hand side above (because there were no common factors to cancel in the above).
> difficulties := normal(rhs("));
  difficulties :=

                 6        5         4         3         2
           e (7 e  + 168 e  + 1365 e  + 4368 e  + 5145 e  + 2520 e + 251)
      ------------------------------------------------------------------------
       8       7        6         5         4         3         2
      e  + 30 e  + 333 e  + 1692 e  + 4023 e  + 4470 e  + 2663 e  + 576 e + 36
This is the equation satisfies.
> op(alpha);
                         2         2
                       _Z  e - _Z e  - 6 _Z e - _Z + 7 e

> problems := numer(normal(subs(_Z=difficulties,")));

                                   8        7    10         5       2        3
  problems := 36 e (1 + 22 e + 45 e  - 760 e  + e   - 3132 e  + 45 e  - 760 e

               6       9         4
       + 2258 e  + 22 e  + 2258 e )
Now this expression will be zero precisely when we must re-do the calculation.
> factor(problems);
                                2            2        6
                         36 e (e  + 14 e + 1)  (e - 1)
We see that there are problems if e=1 (which we have already handled) and further difficulties if e is a root of , or .
> e := RootOf(e^2+14*e+1,e);
                                        2
                          e := RootOf(_Z  + 14 _Z + 1)
Now we re-factor the characteristic polynomial, using Maple's evala facility for working with algebraic numbers like e above.
> evala(Factor(chrp));

                                             3             4
                    (7 + lambda) (lambda - 1)  (lambda + 1)
We see that the roots now are simple: -1 (now 4 times instead of 3), -7, and +1. It turns out that for there are no new difficulties, and the eigenvectors can be found as before. However, for , we find that there are only 3 independent eigenvectors: the geometric multiplicity turns out to be only 3. This is a problem well worth uncovering.
> Ap := map(simplify@eval,op(A)):

> R3a := RowEchelon(subs(lambda=-1,op(Ap)),'D3a');

            [                                     2                       ]
            [ 1  0  0  0  0  0   0       RootOf(_Z  + 14 _Z + 1) + 14     ]
            [                                                             ]
            [ 0  1  0  0  0  0   -1                   0                   ]
            [                                                             ]
            [ 0  0  1  0  0  -1  0                    0                   ]
            [                                                             ]
            [                                       2                     ]
            [ 0  0  0  1  0  1   1   - 1/2 RootOf(_Z  + 14 _Z + 1) - 13/2 ]
     R3a := [                                                             ]
            [                                       2                     ]
            [ 0  0  0  0  1  1   1   - 1/2 RootOf(_Z  + 14 _Z + 1) - 13/2 ]
            [                                                             ]
            [ 0  0  0  0  0  0   0                    0                   ]
            [                                                             ]
            [ 0  0  0  0  0  0   0                    0                   ]
            [                                                             ]
            [ 0  0  0  0  0  0   0                    0                   ]
We could have used SPM_quotalias" to clean up that output a bit, as we did earlier for SPM_quotalpha". However, this is not too bad.

The rank is 5, so there are only 3 free parameters for finding the eigenvectors. In fact, we could have predicted this, because our earlier solution for got all the eigenvectors that there could be.



next up previous
Next: Remarks on the Up: Eigenvalue Problems. Previous: Symbolic Computation With



Robert Corless
Wed Jan 31 11:33:59 EST 1996