Doc/Articles/Play141
At Play With J ... Table of Contents ... Previous Chapter ... Next Chapter
15. Oh, No, Not Eigenvalues Again!
. By Eugene McDonnell. First published in Vector, 14, 1, (July 1997), 135-139.
I can't explain why it is that I keep running into problems associated with eigenvalues. I don't seek them out, and have no interest in them, but there it is--they keep cropping up before me and I have somehow to find a way to drive a stake through their hearts before I can go on to something else.
When the eigenvalue problem last loomed before me, I found that I had to go back to basics in order to come to terms with it. In the course of doing so, I happened upon a technique in J for finding eigenvalues which is eminently satisfying pedagogically, since it can be used for small matrices to show all the theory of eigenvalues, even though it is staggeringly inefficient, becoming dreadfully slow for matrices of size 5 or 6 or larger. It is its pedagogical merit that I commend to you.
The basics are simple. If A is a square numeric matrix, we replace each diagonal atom aii with the two-atom list aii _1. What we are doing is replacing the problem of evaluating the determinant of a numeric matrix with that of evaluating the determinant of a matrix of polynomials. Of course, the non-diagonal terms are the simplest of polynomials, that is, constants, but the diagonal terms are all linear polynomials. Our problem is to evaluate the determinant of this polynomial matrix. In J we evaluate the determinant of a numeric matrix with the monad of the dot conjunction:
det =: -/ . * NB. (1)
For example,
]m =: 2 3,:5 8 2 3 5 8 det m 1
The verb det is not directly suited to our purpose, but with a few maneuvers we'll get where we want to go. First, let's take a suitable matrix for our example (the example and much of the approach I have borrowed from Cornelius Lanczos's book, Applied Analysis, Prentice-Hall, 1956, chapter II):
]A =: 33 16 72 , _24 _10 _57 ,: _8 _4 _17 33 16 72 _24 _10 _57 _8 _4 _17
We can obtain the diagonal atoms of A using the monad d [JP 3.B.m4, meaning phrase m4 from section B of chapter 3 of the book J Phrases].
d =: (<0 1)&|: NB. JP 3.B.m4 ]a0 =: d A 33 _10 _17
We append _1 to each of these using J's stitch (,.):
]a1 =: a0 ,. _1 33 _1 _10 _1 _17 _1
In order to be able to work with a matrix in which a single number is replaced by a list of two numbers, it will be necessary to box the rows of the 2-column matrix above, using monad B1 [JP 1.C.m11]
B1 =: <"1 NB. JP 1.C.m11 ]a2 =: B1 a1 +-----+------+------+ |33 _1|_10 _1|_17 _1| +-----+------+------+
In order to amend the diagonal atoms of A, we'll box the atoms of matrix A first, using monad B0 [JP 1.C.m10]:
B0 =: <"0 NB. JP 1.C.m10 ]a3 =: B0 A +---+---+---+ |33 |16 |72 | +---+---+---+ |_24|_10|_57| +---+---+---+ |_8 |_4 |_17| +---+---+---+
The amend adverb in J needs the indices of the places to be amended. We obtain these using the monad d again, modified by the adverb IR [JP 3.B.a3].
IR=: @(i.@$@]) NB. JP 3.B.a3 ]a4 =: (d IR) A 0 4 8
We can use amend now to obtain the matrix we want:
]a5 =: a2 d IR } a3 +-----+------+------+ |33 _1|16 |72 | +-----+------+------+ |_24 |_10 _1|_57 | +-----+------+------+ |_8 |_4 |_17 _1| +-----+------+------+
For use later, we'll collect the steps taken so far to form the verb db:
db =: B1@(_1 ,.~ d) d IR } B0
Now we are in position to evaluate the determinant of this matrix. We can't use the det verb from above, since it can only deal with matrices whose atoms are simple numbers. My first thought was to use the definition of the Determinant conjunction u . v from the J Introduction and Dictionary:
u =: -/
v =: *
DET =: v/@,`({."1 u .v$:@minors)@.(1&<@{:@$)"2
minors =: }."1@(1&([\.))
replacing the definitions of u and v by polynomial difference and polynomial product, dyads dif and ppr [JP 9.C.d1 and 9.C.d2]. However, I hadn't proceeded far when I suddenly realized that, since DET was the definition of the monadic form of the dot conjunction, I ought to be able to use dot itself directly, and not mess with its sybilline definition. The first thing I did was to provide myself with modified forms of dif and ppr, ones which operated on boxed atoms:
dif =: -/@,: NB. JP 9.C.d1 ppr =: +//.@(*/) NB. JP 9.C.d2 difb =: dif&.> pprb =: ppr&.>
Taking this tack, I could now write, in direct analogy with (1),
detp =: difb/ . pprb
and then apply this to matrix a5:
detp a5 +----------+ |6 _11 6 _1| +----------+
I wanted an open, not a boxed result, so modified the definition of detp:
detp =: >@(difb/ . pprb) ]a6 =: detp a5 6 _11 6 _1
This looked promising, but was it indeed the characteristic polynomial of my matrix? One way to find out was to use the Cayley-Hamilton theorem, which says that a matrix satisfies its own characteristic equation. To see whether this was so, I found the first four powers of A, using the monad I [JP 9.A.m0] to form the 0th power of the matrix, that is, the identity matrix:
I =: =@i.@# NB. JP 9.A.m0 ip =: +/ . * NB. inner product ]a7 =: A&ip^:(i.@>:@#`]) I A 1 0 0 0 1 0 0 0 1 33 16 72 _24 _10 _57 _8 _4 _17 129 80 240 _96 _56 _189 _32 _20 _59 417 304 648 _312 _220 _507 _104 _76 _161
then multiplied these by the supposed coefficients of the characteristic equation:
]a8 =: a6 * a7 6 0 0 0 6 0 0 0 6 _363 _176 _792 264 110 627 88 44 187 774 480 1440 _576 _336 _1134 _192 _120 _354 _417 _304 _648 312 220 507 104 76 161
and lastly, summed a8, trusting that the result would be the zero matrix:
+/a8 0 0 0 0 0 0 0 0 0
From here it was a short step to the eigenvalues: they are the roots of this polynomial, so I could use J's polynomial primitive (p.):
p. detp a5 +--+-----+ |_1|3 2 1| +--+-----+
The result I wanted was the open (>) of the tail ({:) of this one, so I made another modification:
]a9 =: >@{:@p.@detp a5
3 2 1
The parts can now be assembled to give the eigenvalue finder cm:
cm =: > @ {: @ p. @ detp @ db
cm A
3 2 1
If these are the eigenvalues, their product should be equal to the determinant of the matrix.
*/ cm A 6 det A 6
So far, so good. Furthermore, if these are the eigenvalues, if any one of them is subtracted from the diagonals of A, the determinant of the result should be zero. Is this the case?
det A - 3 * I A 0 det A - 2 * I A 0 det A - 1 * I A 0
Yes it is.
The verb cm is a model of the monad of the c. primitive described in the J Introduction and Dictionary 1996 Edition. It differs in that the roots are given in order of descending magnitude, which is how the polynomial (p.) primitive provides them, rather than the ascending order prescribed in the Dictionary. Since c. has not yet been implemented, it's anyone's guess how p. and c. will be reconciled. I've brought this matter to the attention of the proper authorities, so they do at least know that the problem exists.
Vector Vol.14 No.1