65 RAČUNALNIŠKI PODPROGRAMI INFORMATICA 2/91 ZA RAČUNANJE LASTNIH VREDNOSTI IN LASTNIH VEKTORJEV Keywords: principal components analysis, Vesna Čančer eigenvalues, eigenvectors, diagonalization Ekonomsko-poslovna fakulteta v Mariboru method, subroutine Razlagova 14, 62000 Maribor V članku obravnavamo računalniške podprograme za računanje lastnih vrednosti in lastnih vektorjev simetrične korelacljske matrike, ki Jih potrebujemo pri analizi glavnih komponent. V nJem Je uporabljena modifikacija Jacobi jevega postopka, imenovana "Pope in Tompkins", ki Je primerna za računalniško programiranje. Čeprav se uporabljajo hitrejše metode s tridlagonalizirano korelacljsko matriko, se Jacobijeva metoda Se pojavlja v sodobnih podprogramsklh paketih, saj omogoča natančne Izračune lastnih vrednosti in nJim pripadajočih lastnih vektorjev. Delovanje podprograma smo ponazorili s primerom in izpisali končni rezultat. ABSTRACT The following paper deals with the fortran computer subroutines for computing eigenvalues and eigenvectors of a real symmetric matrix, which are necessary for the principal components analysis. Because of its convenience for computer programming, diagonalization method originated by Jacobi and modified by the iterative "Pope and Tompkins" technique is used. In spite of some quicker methods, based on the threediagonallzation matrix, Jacobi's method still appears in the modern subroutine packages because it provides precise definitions of eigenvalues and associated eigenvectors. A numeric example illustrating the application of the described procedure and given computer subroutine is presented. 1. UVOD Analiza glavnih komponent Je statistična tehnika, s katero n statističnih znakov x , x , 1 2 .... linearno in ortogonalno nadomestimo z enakim številom nekoreliranlh glavnih komponent yt , y2..........[3], [8], Transformacija temelji na iskanju lastnih vrednosti in lastnih vektorjev kovariariCne matrike ali, kot v našem primeru, korelacljske matrike A, ki je kovarianCna matrika standardiziranih spremenljivk x , x , ..., x 13J. 12 rt [4], [6], Vzemimo n-razse2ni vektor z elementi r , r , i s . . . , r in poskušajmo rešiti matrično enačbo 1 n 2 n ct.o kjer je A matrika korelacijskih koeficientov In K poljubno Število [8]. To enačbo Je mogoče zapisati v obliki sistema n linearnih enačb, ki Jih uredimo tako, da vse Člene z neznankami r^, r^, ..., r prenesemo na levo stran. Dobimo homogen sistem enačb; desne strani vseh enačb so namreč enake nic [8], [10J. Tak sistem ima pri poljubni vrednosti \ trivialno rešitev r - r =..."» r - 0. Netrivialna 12 r, , rešitev . eksistlra kvečjemu tedaj, kadar je 66 determinanta koeficientov A •» 0 vzamemo tako vrednost, da enačbe sistem, dobimo neskončno mnogo [10]. Ce za X tvorijo odvisen rešitev. Sistem homogenih enačb Je torej determinanta sistema enaka nlc: odvisen, ce je l-x 1 2 1-X 1-X (1.2> To Je karakteristična enačba korelacijske matrike, ki Jo lahko zapišemo v obliki algebraJske enačbe n-te stopnje [8], Rešitve karakteristične enačbe <1.25 K^, K^, .... imenujemo lastne vrednosti korelacijske matrike. Predpostavimo, da so lastne vrednosti realne in različne. Uredimo jih po velikosti: X X > ... > X . 2 n <1.3> Vsaki lastni vrednosti X , i ■» 1, 2, ..., n, v pripada lastni vektor r« , ki Je rešitev matrične enačbe: Ar »X r . v v v. Lastnosti lastnih vrednosti in lastnih vektorjev korelacijske matrike so podrobneje opisane npr. v [3]. [4], [8] in [10]. S podprogramom EIGEN [5] se računajo lastne vrednosti in lastni vektorji simetrične korelacijske matrike, katere elementi so v podprogramu urejeni v vektor, označen s formalno spremenljivko A. V tej spremenljivki se izračunajo tudi lastne vrednosti. Tvori se matrika Lastnih vektorjev, katere elementi so v podprogramu urejeni v vektorju R. 2. OSNOVE UPORABLJENE METODE Uporabljena je za računalniško programiranje primerna modifikacija JacobiJeve metode za določanje lastnih vrednosti, t. J. algoritem Pope in Tompkins [4], ki ga Je za računalnike priredil von Neumann. Bistvo JacobiJevega postopka za ugotavljanje lastnih vrednosti Je, da se z zaporedjem ortogonalnih transformad J korelacijska matrika pretvori v diagonalno matriko lastnih vrednosti: a R' A R l l = R' A R 2 12 A ■= R'A R => A, g 9 g-' g kjer Je A simetrična korelacijska matrika, R transformacijska matrika ortogonalnih transformacij in A diagonalna matrika lastnih vrednosti. Z vsako ortogonalno transformacijo lahko en nedlagonalni element reduciramo na nic. Proces se prične z redukcijo največjega nediagonalnega elementa, nato se reducira drugi največji nedlagonalni element itd., dokler niso vsi nedlagonalni elementi reducirani na nlc. V vsakem zaporedju moramo poznati elementarno transformacijsko matriko, s katero bomo izbrani element v matriki reducirali na nic [10], [11]: cos © -sin G 0 Sine O cos 0 0 0 1 <2.1> Kot 0, za katerega moramo zarotirati matriko A, da bo element alm postal nic, Je z izrazom: korelaci jsko podan tang 2 © = -2 a, / sin © cos © + a. Ccos © - sin ©> »0. mm lm <2.3) Upoštevajmo vrednosti za sin 2 0 in cos 2 0: <2.4) 1/2 sin 2 © + a cos 2 © 11 mm lm iz cesar sledi: 0, -2 a. ' Nato iz izraza za t.angens dobimo vrednosti sin © in cos ©, t.J. elemente transformacijske matrike R. Bistvena razlika med izvirno JacobiJevo metodo in algoritmom Pope in Tompkins Je v tem, da se z njim na nic ne reducira le en element, pac pa se z določitvijo minimalne meje tolerance zaporedno eliminirajo vsi nediagonalni elementi, ki so večji od te meje. Nato se določi nova minimalna meja tolerance, ki Je nižja od prejšnje meje. Proces se nadaljuje vse do konCne meje tolerance napake za lastne vrednosti. Algoritem tvori osem korakov: V 1. koraku poiščemo vsote kvadratov vseh nediagonalnih elementov korelacijske matrike in kvadratni koren te vsote, t. J. začetno normo E 2 a mm <2.8> C2.9) C2.10) sin © <2. 11 > -/(2C1 + -/ci - ■ w*)) cos © - /o - sin2©) C2. 12) B ■> a. l cos 0 - a. v m sin © C2. 13) C = a i l sin © + a, L rn cos © <2 . 14) B = r i l cos © - r l m sin © (2. . iS) r i m = r sin i l © + r. L T cos © n <2. . 16) i l \l B <2.17) a cos2© + a sin2© - 2a. sin © cos © l l mm Im <2.18> mm l l sinZ© + a cos2© + 2a, sin © cos © a ) sin 0 cos , + mm + a Ccos © - sin ©> <2.19) <2.20) 3. OPIS DELOVANJA PODPROGRAMA EIGEN S PRIMEROM V glavnem programu [2] kličemo podprogram EIGEN z ukazom CALL, s katerim se dejanski parametri iz glavnega programa priredijo formalnim parametrom v podprogramu: formalni pomen Saratnel n poaprog. A vhodna simetrična korelaci jska matri- ka, urejena kot vektor, ki se med računanjem v podprogramu uniči; izhodni vektor lastnih vrednosti, ki so elementi matrike lastnih vrednosti R izhodna matrika lastnih vektorjev, urejena kot vektor N rang matrik A in R, enak redu teh ma- trik zaradi neodvisnosti statističnih znakov MV vhodna spremenljivka za EIGEN, ki Ji priredimo vrednost z ukazom DATA v glavnem programu in Je lahko 0 za računanje lastnih vrednosti in ve k t o r Je v ali _1 za računanje zgolj lastnih vrednosti V glavnem programu smo spremenljivki MV priredili vrednost 0. V EIGEN-u niso potrebni drugI podprogrami ali funkcijski podprogrami. Za ponazoritev delovanja podprograma bomo navedli. vmesne in izpisali končne rezultate. Za lažje razumevanje uporabljenega algoritma izplslmo podprogram EIGEN [S]. SUBROUTINE EIGEN(A,R,N,MV) DIMENSION A( 1),R(1) C TVORJENJE ENOTSKE MATfelKE 5 RANGE=1.OE-6 IF(MV-l) 10,25,10 10 IQ=-N DO 20 J = 1, N IQ=IQ+N DO 20 1=1,N IJ=IQ+I R(IJ )=0.0 IF(I-J) 20,15,20 15 R(IJ ) = 1.0 20 C0NTINUE C RAČUNANJE ZAČETNE IN KONČNE NORME (ANORM IN C ANRMX) 25 ANORM=0.0 DO 35 I-l.N DO 35 J=I,N IF(I-J) 30,35,30 30 IA=I+(J#J-J)/2 ANORM=ANORM+A(IA)*A(IA) 35 C0NTINUE IF(ANORM) 165,165,40 40 ANORM=1.414*SQRT(ANORM) ANRHX=AN0RM*RANGE/FL0AT(N) C INICIALIZACIJA INDIKATORJA IN RAČUNANJE C INPUTA (THR) IND = 0 68 45 50 55 60 62 65 68 70 75 78 80 85 90 95 100 105 110 115 120 125 130 135 140 145 150 155 160 165 170 175 THR=AN0RM THR=THR/FL0AT(N) L = 1 M=L+1 RAČUNANJE SIN IN COS MQ=(M*M-M)/2 LQ=(L*L-L)/2 LM=L+MQ IF( ABS(A(LM))-THR) 130,65,65 IND = 1 LL=L+LQ MM=M+MQ X=0.5*(A(LL)-A(MM>) Y=-A(LM)/ SQRT(A(LH)*A(LM)+X*X) IF))) SINX2=SINX*SINX COSX= SQRT(1.0-SINX2) COSX2=COSX*COSX SINCS =SINX#C0SX ROTIRANJE L IN M STOLPCEV ILQ=N*(L-1) IMQ=N#(M-1) DO 125 1 = 1,N IQ=(1*1-1)/2 IF(I-L) 80,115,80 IF(I-M) 85,115,90 IM=I+MQ GO TO 95 IM=M+IQ IF(I-L) 100,105,105 IL=I+LQ GO TO 110 IL=L+IQ X=A*SINX2+A(MM)*C0SX2+X A(LM)=(A se v podprogramu EIOEN priredi kot vektor v formalni spremenljivki A z naslednjim zaporedjem elementov: 1., 0.887, 1., -0.125, 0.307, 1., torej Je njihovo Število NOCN+l)/^ « 6. V 1. zapisu Je definiran' podprogram EIGEN in izhodni spremenljivki A in R ter vhodne spremenljivke A, N ln MV. V podprogramu Je upoštevana enojna natančnost. Ce bi želeli dvojno natančnost, bi to definirali z ukazom DOUBLE PRECISION A,R,ANORM,ANRMX,THR,X,Y,SINX,SIX2, COSX,COSX2,SINCS ,RANQE ,DSQRT,D ABS Tudi f ortranske funkcije morajo biti dvojne natančnosti. SQRT v ukazih 40,68,75 ln 78 moramo zamenjati z DSQRT. ABS v ukazu 62 moramo zamenjati v D ABS. Popraviti moramo tudi konstanto v ukazu 5, in sicer na 1.0D-12. Z ukazi 10-1 do 20 se v primeru, ko smo v glavnem programu spremenljivki MV dodelili vrednost O, tvori enotska transformacijska matrika <2.1), katere elementi so urejeni v vektor RCIJ). Izvede se 3. korak algoritma Pope in Tompkins. Z ukazi 25 do 40 se tvori zafietna norma ANORM c.2.6). Izvede se 1. korak algoritma Pope in Tompkins. V našem primeru je začetna norma 1.3389340. V ukazu 40+1 se tvori končna meja tolerance napake za lastne vrednosti oz. konCna norma ANRMX C2.7). Izvaja se torej 2. korak obravnavanega algoritma. V naSem primeru je končna norma U.00U0004. Z ukazi 45-2 do 55 poteka lnicializaclja Indikatorjev IND in določanje začetne minimalne meje tolerance, kar sodi v 4. korak v algoritmu Pope ln Tompkins. V našem primeru se je izračunalo 14 minimalnih mej tolerance THR; prva Je bila 0.4463114, zadnja pa 0.0000003. Z ukazi od 60 do 78+2 se računajo sinusi ln kosinusi, nato pa se do 130-1 zamenjujejo L in M stolpci. V vsaki iteraclji se . namreč 'V') 69 selekcionirajo nediagonalni elementi. To sodi v S., 6. in 7. korak obravnavanega algoritma. Iz ukaza 62 Je razvidno: Ce je nediagonalni element, korelacijske matrike po absolutni vrednosti manjši od prej določene minimalne meje tolerance THR, ni reduciran na nic; Ce so nediagonalni elementi po absolutni vrednosti enaki ali večji od prej določene minimalne meje tolerance, pa se začne reduciranje tega elementa na nic. V našem primeru Je v prvi iteraciji A<2), ki Je 0.887, večji od THR, ki je 0.446, zato sledi reduciranje A<2) na nic. V ukazu 68-1 se izračuna X po <2.95. X predstavlja del Izraza <2.4). V ukazu 68 se po <2.105, v katerem nastopa <2.8), izračuna Y, ki je potreben pri računanju vrednosti elementov transformacijske matrike R in A: sin © in cos ©. V ukazu 75 se po <2.11) izračuna element transformacijske matrike R sin © in se zapomni v formalni spremenljivki SXNX. V naslednjem ukazu se izračuna sin2© in se zapomni v formalni spremenljivki SINX2. V ukazu 78 se po <2.12) izračuna cos © in se zapomni v COSX. V naslednjih dveh ukazih se izračunata Se cos2© in sinScosO in se zapomnita v spremenljivkah COSX2 in SINCS. V 6. koraku se določijo elementi zaporednih transformacijskih matrik po vsaki iteraciji. V ukazu 110 se izračuna X po <2.13). V naslednjem ukazu se izračuna A po <2.14). Nato se A izračuna R in nato po <2.17) Se R