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) ------------------------------------------------------------------------------- eIf

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

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]) eNow check for special cases -- we look at the determinant

> normal(D1); 4Nonzero, independent of

> 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 ------------ eZero 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`

`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_quot`

RootOf" and do just one eigenvector calculation.
Note that `SPM_quot`

alpha" 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_quot`

normal" 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_quot`

evala@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_quot`

simplify", which knows about RootOf and can simplify things like
and to expressions linear in , and `SPM_quot`

normal", 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

> 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 eNow 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_quot`

solve" 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 + 36is a function of

`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 + 36This 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 := 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 > evala(Factor(chrp)); 3 4 (7 + lambda) (lambda - 1) (lambda + 1)We see that the roots now are simple:

> 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_quot`

alias" to clean up that output
a bit, as we did earlier for `SPM_quot`

alpha". 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.

Wed Jan 31 11:33:59 EST 1996