%%HP: T(3)A(R)F(.); @ Author: Ted A Smith @ Date: Sun July 21, 1991 @ @ Here is my second cut at EigenValue/EigenVector decomposition @ which is not restricted to real symmetric input matricies. These @ routines have been optimized (and accuracy improved) compared to my @ first try dated Jul 18 1991. @ @ Checksum of this directory on the stack: #2B90h, size 2982.5 @ Checksum with only Eigen, EFn1 and EFn2: #F1BBh, size 1103.5 @ @ The general eigenvalue/eigenvector problem is hard and there are @ many posible pathologies if the input matrix is not real symmetric. @ (For theory and examples see the references below.) @ @ The directory consists of 3 main routines: @ @ Eigen Given M returns V and D, where V are the eigenvectors @ and D are the eigenvalues @ EFn1 Given M and f returns f(M) defined by V*f(D)*V^-1 @ EFn2 Given V, D and f returns V*f(D)*V^-1 @ @ For example, in analogy with SIN(x)^2+COS(x)^2=1: @ @ A Eigen DUP2 \<< SIN \>> EFn2 DUP * @ ROT ROT \<< COS \>> EFn2 DUP * + A IDN - ABS @ @ will return a small number for almost any square input matrix A. @ @ Eigen also displays the matrix after each rotation if user flag 1 is @ set. Two comments in the code indicate which lines to delete if this @ feature is not desired. @ @ There are two test routines: @ @ Tst Calculate ABS(DET(M-D*I)) and ABS(M-EFn1(M, \<< \>>)) @ Tst1 Apply SIN(x)^2+COS(x)^2=1 test to a random matrix. @ @ There are 9 example matricies A, B, C, D, E, F, G, H and J as well as @ the 3x3 identity matrix I. The comments in the code give the true @ eigenvalues for some of these matricies. @ @ Here following is the theory of operation: @ @ The QR algorithm is, in theory, a better method; but it is much @ more complicated. The virtue of this implementation is its simplicity. @ It is based on the Jacobi method, which uses trignometric rotations to @ succesively zero the off-diagonal elements of the input matrix. When @ the off-diagonal elements are close enough to zero the diagonal @ elements approximate the eigenvalues of the original matrix. In @ addition, accumulating the product of the rotations gives the matrix of @ eigenvalues. On real symmetric matricies the Jacobi method is @ idiotproof. @ @ I have extended the algorithm by using hyperbolic "rotations" to @ drive the off-diagonal elements tword each other, thereby allowing the @ trignometric transform to be more effective. I have no proof that this @ converges! I have not been able to find a reference for this method. @ The idea of using hyperbolic transforms was inspired by chapter 17 of @ reference 3, but this implementation bears no other obvious similarity. @ @ To illustrate the method let's apply it to the 2 x 2 matrix: @ @ [[a b] @ [c d]] @ @ The standard Jacobi rotation is: @ @ [[ cos(t) sin(t) ] [[a b] [[ cos(t) -sin(t) ] @ [-sin(t) cos(t) ]] [c d]] [ sin(t) cos(t) ]] @ @ where t=ATAN((b+c)/(a-d))/2. @ @ For real a, c and d and b=c this returns: @ @ / a+d+SQRT(amd^2+4*c^2)*SIGN(amd) \ @ | -------------------------------- 0 | @ | 2 | @ | a+d-SQRT(amd^2+4*c^2)*SIGN(amd) | @ | 0 ------------------------------- | @ \ 2 / @ @ where amd=a-d. @ @ To drive the off-diagonal elements together we will first apply the @ hyperbolic transform: @ @ [[ cosh(t) sinh(t) ] [[a b] [[ cosh(t) -sinh(t) ] @ [ sinh(t) cosh(t) ]] [c d]] [-sinh(t) cosh(t) ]] @ @ Where t=ATANH((b-c)/(a-d))/2. @ @ Once again, for real inputs (with fortunate signs) this returns: @ @ / a+d+SQRT(amd-bmc)*SQRT(amd+bmc) b+c \ @ | ------------------------------- ----- | @ | 2 2 | @ | | @ | b+c a+d-SQRT(amd-bmc)*SQRT(amd+bmc) | @ | ----- ------------------------------- | @ \ 2 2 / @ @ where amd=a-d and bmc=b-c. @ @ The transforms work with complex numbers, however the on-diagonal @ results get much grosser. @ @ To save matrix math the code utilizes one transform which is the @ product of the trig and hyperbolic transforms: @ @ [[ cb*ca+sb*sa cb*sa+ca*sb ] [[a b] [[ cb*ca-sb*sa -cb*sa-ca*sb ] @ [ cb*sa-ca*sb cb*ca-sb*sa ]] [c d]] [ ca*sb-cb*sa cb*ca+sb*sa ]] @ @ where cb=COS(beta), sb=SIN(beta), ca=COSH(alpha) and sa=SINH(alpha) @ and alpha=ATANH(bmc/amd)/2 @ and beta=ATAN((b+c)/(amd*SQRT(1-bmc/amd)*SQRT(1+bmc/amd)))/2 @ and amd=a-d and bmc=b-c. @ @ (The seemingly redundant math assures that signs, etc. are correct for @ complex inputs.) @ @ I had implemented the calculation of the xform matricies without @ using explicit trig and hyperbolic functions, but getting all of the @ edge conditions right is a nightmare (at least for me). Since those @ implementations only sped up the aggragate algorithm by only 2-3% I @ opted for simplicity. @ @ Known limitations: @ Defective matricies cannot be eigenvalue/eigenvector decomposed. @ Hence, this and any other implementations are doomed in the general @ case. (See ref 2 pages 353... for an explanation.) @ @ Although it would be nice, I do not yet perform orthogonalization of @ the resultant eigenvectors, so matricies with repeated eigenvalues @ usualy do not yield orthogonal eigenvectors. @ @ If the input matrix is Real Symmetric, I still use complex math. @ There is a posibility that roundoff errors will cause complex @ numbers to spontainiously errupt causing an error if I were to @ try to stuff them into a real array. Try example matrix H. @ It would be nice to take advantage of real only math to speed @ up the algorithm. @ @ Certain input matricies will contain values which will cause the trig @ xform to attempt to evaluate ATAN(i) or the hyperbolic xform to @ evaluate ATANH(1). This indicates that no possible rotation will @ work in that context. Since no rotation will work, I do no @ rotation. If no other rotation happens before we get back to this @ same position on the next pass, we won't converge. Example matrix @ C exibits this problem during the first pass when rotating at @ position 1, 2. @ @ Suspected problems: @ The test for termination is a hack: I quit when ABS(M) is not much @ different than ABS(diag(M)), it seems to work fine. @ @ @ References: @ @ 1) Matrix Computations, 2nd Edition, 1989 @ Golub, Gene H. and Van Loan, Charles F. @ The John Hopkins University Press @ 701 West 40th Street @ Baltimore, Maryland 21211 @ ISBN 0-8018-3772-3 @ ISBN 0-8018-3772-1 (pbk.) @ @ 2) Numerical Recipes (in C, Fortran or Pascal) @ The Art of Scientific Computing @ Press, William H., Flannery, Brian P., @ Teukolsky, Saul A. and Vetterling, William T. @ Cambridge University Press @ 510 North Avenue @ New Rochelle, NY 10801 @ ISBN 0-521-35465-X @ @ 3) Linear Algebra, vol. II of Handbook for Automatic Computation, 1971 @ Wilkinson, J.H. and Reinsch, C. @ Springer-Verlag @ New York DIR Eigen \<< 1 IF FS? THEN CLLCD END @ Nuke this line for no display RCLF SWAP -55 SF -21 CF -22 SF DUP RE SWAP IM R\->C DUP IDN SWAP DUP SIZE 1 GET \-> d \<< DO 1 d 1 - FOR i i 1 + d FOR j 1 IF FS? THEN DUP 1 DISP END @ Nuke this line for no display DUP i j 2 \->LIST GET OVER j i 2 \->LIST GET DUP2 + ROT ROT - 3 PICK i DUP 2 \->LIST GET 4 PICK j DUP 2 \->LIST GET - SWAP OVER IFERR / THEN 0 END ROT ROT 1 4 PICK DUP2 + \v/ ROT ROT - \v/ * * IFERR / THEN 0 END ATAN 2 / IF DUP ABS 8 > THEN DROP 0 END DUP COS SWAP SIN ROT ATANH 2 / IF DUP ABS 8 > THEN DROP 0 END DUP COSH SWAP SINH 5 PICK IDN DUP \-> T TT \<< 4 DUPN SWAP 4 ROLL * SWAP ROT * DUP2 + DUP 'T(i,i)' STO 'TT(j,j)' STO - DUP 'T(j,j)' STO 'TT(i,i)' STO 4 ROLL * SWAP ROT * DUP2 + DUP 'T(i,j)' STO NEG 'TT(i,j)' STO - DUP 'T(j,i)' STO NEG 'TT(j,i)' STO T ROT * T ROT * TT * \>> NEXT NEXT UNTIL 0 1 d FOR i OVER { i i } GET DUP CONJ * RE + NEXT \v/ OVER ABS %CH ABS .000000001 \<= END 1 d FOR i DUP i DUP 2 \->LIST GET SWAP NEXT DROP d \->ARRY \>> ROT STOF \>> EFn1 @ Apply function to arbitrary matrix \<< SWAP Eigen ROT EFn2 \>> EFn2 @ Apply function f to decomposed matrix \<< \-> v f @ If M=V*D*V^-1, F(M) is V*F(D)*V^-1 \<< DUP 0 CON 1 OVER SIZE 1 GET @ Loop over the elements of D FOR i { i i } v i GET @ Eval f at each element of D f EVAL PUT @ And build F(D) NEXT OVER / SWAP * @ Do V*F(D)*V^-1 \>> \>> Tst @ Check Eigen(M) three ways \<< DUP Eigen DUP SIZE 1 GET \-> M T D d \<< 1 d @ Test Eigenvalues by doing DET(M-v*I) FOR i @ for each eigenvalue v M DUP IDN D i GET * - DET NEXT d \->ARRY ABS "Det" \->TAG T DUP TRN CONJ * DUP IDN - @ Check orthogonality of eigenvectors ABS "Orth" \->TAG @ with ABS of DOT of all pairings M T D \<< \>> EFn2 - ABS @ Check eigenvectors with "Res" \->TAG 3 \->LIST @ ABS(M - V * [v] * INV(V)) \>> \>> Tst1 @ Test by calculating SIN(M)^2+COS(M)^2 \<< @ for a random matrix M. The result RAND RAND RAND RAND RAND @ should be close to the identity RAND RAND RAND RAND @ matrix. { 3 3 } \->ARRY RAND RAND RAND RAND RAND RAND RAND RAND RAND { 3 3 } \->ARRY R\->C DUP Eigen DUP2 \<< SIN \>> EFn2 DUP * ROT ROT \<< COS \>> EFn2 DUP * + \>> I [[ 1 0 0 ] [ 0 1 0 ] [ 0 0 1 ]] A [[ 1 2 3 ] @ 11/3-2*SQRT(133)*SIN(ATAN(433*SQRT(3)/117)/3)/3 [ 2 4 5 ] @ 11/3+2*SQRT(133)*SIN(ATAN(433*SQRT(3)/117)/3+pi/3)/3 [ 3 5 6 ]] @ 11/3-2*SQRT(133)*COS(ATAN(433*SQRT(3)/117)/3+pi/6)/3 B [[ 2 0 3 ] @ Exact eigenvalues are 0, 2 and 5 [ -1 -1 -3 ] [ 1 3 6 ]] C [[ 1 0 .01 ] @ Eigenvalues are 11/10 [ .1 1 0 ] @ 19/20-SQRT(3)*i/20 [ 0 1 1 ]] @ 19/20+SQRT(3)*i/20 D [[ -261 209 -49] @ Eigenvalues are 3, 4 and 10 [ -530 422 -98] [ -800 631 -144]] E [[ (1,1) (0,0) (.01,.01) ] @ Eigenvalues are (11/10)*(1+i) [ (.1,.1) (1,1) (0,0) ] @ (19/20-SQRT(3)*i/20)*(1+i) [ (0,0) (1,1) (1,1) ]] @ (19/20+SQRT(3)*i/20)*(1+i) F [[ (0,5) (1,2) (2,3) (-3,6) (6,0) ] [ (-1,2) (0,6) (4,5) (-3,-2) (5,0) ] [ (-2,3) (-4,5) (0,7) (3,0) (2,0) ] [ (3,6) (3,-2) (-3,0) (0,-5) (2,1) ] [ (-6,0) (-5,0) (-2,0) (-2,1) (0,2) ]] G [[ 1 2 0 0 ] [ 2 3 4 0 ] [ 0 4 5 6 ] [ 0 0 6 7 ]] H [[ -4 1 1 1 ] @ Eigenvalues are 1, 5, 5 and 5 [ 1 -4 1 1 ] [ 1 1 -4 1 ] [ 1 1 1 -4 ]] J [[ 1 -1 0 0 ] [ -1 2 -1 0 ] [ 0 -1 3 -1 ] [ 0 0 -1 4 ]] END @ Ted A Smith @ PO Box 6308 @ Longmont CO 80501 @ H) (303)651-2092 @ W) (303)665-5492 @ Jul 20, 1991 @ akcs.tasmith@@hpcvbbs.hp.cv.com