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.