"Supported in part by the Advanced Research Projects Agency of the Department of Defense and ... monitored by the U.S. Army Research Office-Durham under contract no. DAHCO4-72-C-0001."Bibliography: p. 8

Lermit, J.,

Noh, K.,

Sameh, Ahmed

IDEALS @ Illinois

On the intermediate eigenvalues of symmetric sparse matrices, CAC No. 91

Illinois Digital Environment for Access to Learning and Scholarship Repository

UNIVERSITY OFILLINO'S • IHRARYAT URBANAC: wvlPAlGNENGINEERINGNOTICE: Return or renew all Library Materialsl The Minimum Fee foreach Lost Book is $50.00. 11141 e\The person charging this matei4*l fe rtS||^sible forits return to the library from which it was withdrawnon or before the Latest Date stamped below.!,b^9ks alJP' r^teQns for discipli-' ''ersity.UNIVERSITY ILLINOIS LIBRARY URBANA-CHAMPAIGNNOV i J nnNOV 2 2 neniL161—O-1096;ONFERENCEENGINEERING UBRabvUNlVE«S,Ty OF So'sCenter for Advanced CoriMSSlBSKSlBS^A^iUNIVERSITY OF ILLINOIS AT URBANA-CHAMPAIGNURBANA. ILLINOIS 618012/^—/ ^^uiucDcs-R-73-610CAC Document No. 91ON THE INTERMEDIATE EIGENVALUES OFSYMMETRIC SPARSE MATRICESBy /iA. Sameh, J. Lermlt, and K. NohOctober 1973(Revised February 197^)Digitized by the Internet Archivein 2012 with funding fromUniversity of Illinois Urbana-Champaignhttp://archive.org/details/onintermediateei91sameUIUCDCS-R-T3-610CAC Document No. 91ON THE INTERMEDIATE EIGENVALUES OFSYMMETRIC SPARSE MATRICESByA. Sameh, J. Lermit, K. NohCenter for Advanced ComputationUniversity of Illinois at Urbana-ChampaignUrbana, Illinois 618OIFebruary, 19T^(Revised)This work was supported in part by the Advanced Research ProjectsAgency of the Department of Defense and was monitored by the U.S.Army Research Office-Durham under Contract No. DAHC0l|-T2-C-0001.ABSTRACTAn algorithm has been developed for finding the eigenvalues of asymmetric matrix A in a given interval [a, b] and the correspondingeigenvectors using a modification of the method of simultaneous itera-tions with the same favorable convergence properties. The techniqueis most suitable for large sparse matrices and can be effectively im-plemented on a parallel computer such as the ILLIAC IV.ACKNOWLEDGEMENTWe thank Professor Geneva Belford for her constructive commentsand discussions. We also thank Lois Pelczar for typing the manu-script.IntroductionConsider the eigenvalue protlem Au = Xu where A is a real symmetricmatrix of order n. If n is large and A is sparse, then an iterativemethod for finding the eigenvalues and eigenvectors which uses A only asan operator may be superior to a method which reduces A to a condensedform (e.g., tridiagonal) using orthogonal transformations resulting in afar denser matrix exceeding storage capacity. This superiority is quitedemonstrable if only a few of the eigenvalues and the correspondingeigenvectors are required. Bauer's method of Simultaneous Iterations[l, 2] for finding a few of the leading eigenvalues and eigenvectors isone of such methods. In this paper we modify his method to find theeigenvalues of A in some interval [a, b]. Without loss of generality weassume that A is positive-definite. The problem is therefore to find theq eigenvalues X,Ajt,...,X, ^, wherep p+1 p+q-1c„>X >A^> . . . >X >b>X >X ^ > ... >X , >a>X . > ... >X >c >02— 1— 2— — p-1— p— p+1— — p+q-l — p+q— — n— 1and the corresponding eigenvectors.The AlgorithmSimilar to the method of "Simultaneous Iterations" [2], we proposethe following iterative scheme:(i) Choose an nxq matrix X such thatoX^ X = Io o q^ii) \+k = V^)\ k= 0, 1, 2, ... (1)where B is a polynomial of A to be defined later, and T (B) is themChebyshev polynomial of B of degree m.(iii) Using Gram-Schmidt orthonormalization process we decompose X asJwhere U.U. = I and R. is an upper triancxilar matrix of order q.(iv) LetZ^=AU^. (3)Then form the positive-definite qxq matrixand solve for its eigenvectors Q .J(v) Thus,Go back to (ii) and so on until X AX approaches a diagonal matrix whoseelements are the eigenvalues of A in [a, b]; this occurs when Q ap-proaches the identity matrix.The matrix B in (l) may be obtained as follows:LetA = [2A - (a+b)l]/(b-a); (6)then an eigenvalue of A is given byX = [2A - (a+b)]/(b-a) (?)and those eigenvalues of A in [-1, l] correspond to the eigenvalues ofA in [a, b]. The interval [-c, c] that contains all the eigenvalues ofA is given by,c = max{d , d } (8)where^1 " l^^i " (a+b)|/(b-a)d^ = |2c2 - (a+b)|/(b-a)If some mapping y = f(x) can be found such that f{\) is outside the inter-val [-1, l] for X in [-1, l] and |f(X)|<_l for X outside [-1, l], then thealgorithm described above (i)-(v) will yield the eigenvalues of A in(9)[a, b] and the corresponding eigenvectors.Let us therefore map the interval [-1, l] on the subintervals [-c, -l]and [l, c], one such mapping may be taken as,y = ±/ax+3or X = (y^-6)/a. (lO)X = -1 corresponds to y = ±1, thus1 = -a+3;and X = +1 corresponds to y = ±c, thus2c = a+3.Thereforeand(11)a = |(c2-l) •3 = |(c2+l)Hence the matrix B is taken asB = [2A^ - (c^+l)l]/(c^-l) (12)y = A(B) = [2X^ - (c^+l)]/(c^-l) (l3)i.e., -(^) < y < -1.aWe now introduce the following theorems.Theorem 1 :Let E be the linear space spanning the columns of X . In case ofstable convergence (no reordering of the eigenvalues if the LR-Cholesky de-composition is applied to G ), the angle 4). between the i-th eigenvectoru. and the linear space E = {x|x=T (B)y, yeE }, spanning the columns of-L fj 111 (J ""-LX 5 is asymptotically for j-x» of order 0(q, ) in whichJ -^q. = max { |T (y )|/|T (y )|} (ll;)1 k?fP "^k I I m 1 IieP, where P = {p, p+1, ..., p+q-l}. The proof is quite similarto that of Theorem 2 in [2] and hence will be omitted here (see AppendixI).Theorem 2:The columns of the matrices X as generated by (i)-(v) are such thatl|u. - x.^-^^ll = O(q.J) (15)where q. is as given by (lU). The proof is again similar to that ofTheorem (3) in [2] (see Appendix l).A proper order of the Chebyshev polynomial, m, can be obtained, asin [2], by stipulating that the parallelization of the columns of X,k+mshould not go further than that at most one decimal digit is cancelledout when these columns are orthonormalized, i.e.,l\-i '^^>l '^°'where|y^| = max |y^|.(l6)iePAn approximation of y is obtained from (13) by replacing A by1[2A (G ) - (a+b)]/(b-a) where A (G ) is the maximum eigenvalue of G .^ J ^ J JThus,cosh [(m-l) cosh IPpI] < 10 > andT ^cosh" 10 „ 3 /,„>,m-l < — = — (ITJcosh |y^| cosh" |y^|As soon as one of the eigenvalues of G stagnates, the correspondingJeigenvalue of A can be found within computer accuracy. Once this happenswe can test for acceptance of the corresponding eigenvector. The accep-tance test is that of [3] except that the discoTinting rule is given byf. := f. X /p^tP (18)1 1 T (a. )m 1where a. is an approximation to y, obtained as explained above, and h isthe number of eigenvalues already accepted.Some Numerical Results .The algorithm described above has been implemented in Fortran onUCLA's IBM 360/91 computer. The flow chart is slightly different from thatof Rutishauser [3],zl:=z2:=ks:=0m:=2Perform m steps withChebyshev iterationsks: =ks+mz2:=z2-mYeszl:=zl+lz2:=Uxzlm: =2xniTRED2, IMTQL2 [h]ks:=ks+liAcceptance tests, reduce m if neededTo demonstrate the numerical behavior of the algorithm we tested threesymmetric matrices A , A , and A (shown in Appendix II). In all threeexamples we obtained the required eigenvalues correct to lU decimal placesand the eigenvectors correct to at least 7 decimal places.Example 1The Gerschgorin disks of the Sh'xSk matrix A show that we have 16eigenvalues in the interval (l.i+, 3.6), 8 eigenvalues in the interval(5.^, ^'6)i and i+0 eigenvalues in the interval (8.^, 16.O). We seek toobtain the 8 intermediate eigenvalues and eigenvectors. Two differentintervals [a, b] have been tested [U.O, 8.0], and [5-3, 6.7] with^o"^^kV %2' •••' %8^-TABLE I[a, b] M zl ks[i+.O, 8.0] l6 T lli+[5.3, 6.7] 32 15 U88M: maximum degree of Chebyshevpolynomialszl: niomber of QL stepsks: nximber of iteration stepsTable I shows that, for the same acceptance test, the number of itera-tion steps required (ks ) is smaller when a and b are as far as possiblefrom the eigenvalues to be evaluated (X , A ,_, .... A .), which isp p+1 P+q-1clear from relation (l^) and the preceeding interval transformation.If we choose X = [e,^, e, , ..., e, ] , however, and [a, b]E[5.3, 6.7]the process converges to the required 8 eigenvalues and eigenvectors inonly 265 iteration steps (ks = 265), zl = 10, but M = 58.Example 2The 6'x6 positive-definite matrix A has the simple roots 1, 5, 25,and a triple root at 15. If we seek the triple root 15 and take the in-terval [a, b] as [7, 2k] we require 62 iteration steps, zl = 5, and themaximum degree M of the Chebyshev polynomials is 6. Here we only obtain aninvariant subspace that corresponds to X(A ) = 15.If, however, we seek the four eigenvalues in the interval [2, 2k]we require only 20 iteration steps with the degree of Chebyshev polynomialsnot exceeding 2, and zl = 3. Again we obtain an invariant subspace corre-spending to X(A ) = 15 and the correct eigenvector corresponding to X(A ) = 5.Example 3i+r kTTThe positive-definite matrix A has the eigenvalues X = l6 sin [ / .- \ ],k = 1, 2, ..., n, where n = 6h. There are 26 eigenvalues in (0, 2), 6 in(2, h)y and 32 in (U, l6). Here we woiold like to evaluate the eigenvaluesin the interval [2, h] and the corresponding eigenvectors. Table II showsthe effect of our ass\unption regarding the distribution of eigenvalues onthe number of iteration steps (ks ) required for convergence to the eigen-values and eigenvectors in the interval [2, i+], and degrees of the Chebyshevpolynomials.TABLE IIXoM zl ks(i) [633, ..•' ^38^ 52 9 290(ii) [e^-|_> •••- ^0^ 52 8 237(iii) [e^^, ..••' "37^ 52 21 1312(iv) [e^^, ..'•' ^28^ 52 10 3U3Experiment (iii) indicates that using fewer columns in X than the actualonumber of eigenvalues in [2, h] led to a large ntimber of iteration stepsks = 1312, and zl = 21 to converge to the assumed 3 eigenvalues. Whileexperiment (iv) shows that if we use more columns in X than the actualnumber of eigenvalues, even if we completely misjudge the distribution ofthe eigenvalues in the intervals (O, 2), (2, k), (k^ 16), the n\amber ofiteration steps (ks) and the number of QL steps (zl) required for conver-gence to the 6 eigenvalues and eigenvectors in [2, h] are only slightlymore than those in experiment (i) where we used the true distribution ofthe eigenvalues. Since the nvimber of iterations required appears to belargely independent of the choice of coliomns of the initial matrix X, arandom number generator could be used to generate these.We notice that the above algorithm is quite suitable for parallelcomputations since the major operation here is multiplying a matrix by avector, which can be handled rather efficiently on a parallel computersuch as the ILLIAC IV.REFERENCES[l] Bauer, F. L. , "Das Verfahren der Treppeniteration \ind verwandteVerfahren zur'Losung algebraischer Eigenwertprobleme,"ZAMP, 8 (1957), pp. 21I1-235.[2] Rutishauser, Heinz, "Computational Aspects of F. L. Bauer'sSimultaneoiis Iteration Method," Num. Math., 13 (1969),pp. i+-13.[3] » "Simultaneous Iteration Method for Symmetric Matrices,"Num. Math., 16 (19T0), pp. 205-223.[k] Wilkinson, J. H. , and Reinsch, C, Handbook for Automatic Computation,Vol. II, Linear Algebra, Contributions II/2, II/3. New York:Springer-Verlag, 1971.Appendix IProof of Theorem 1The iterations (i)-(v) are orthogonally invariant, i.e., replacingA by h''^AH = A (where h''^H = I, A = diag (A.)) and X by h"'^X has the ef-•^ 1 o ofeet that all X are replaced by H X , while the G and X are notJ J J Jchanged. Therefore we can ass\ime, without loss of generality, thatA = diag (A^, A^, ..., A^).In case of stable convergence E can be spanned by q vectors^1p-1,1n,l^2'"l,.p-1,2 p-l,qn,2'••1^p+q,l .p+q,2 .p+q,qn,qrow p(A.l)row p+q-110According to (l), E. is spanned byJ\^^^l) -11m p-1 p-1,1m pT "^(y_^ ) X ^ nm.p+q P+q,lT "^(y ) X .m n n,l\^[^1^ \2 \^^) ^,,Tj(y ) x^,_m p-1 p-l>2m p+1JT "^(y J X ,m p-1 p-l,qm p+q-1T "^ (y ) X ... T "^ (y ) xm. P+q P+q,2m. p+q p+q,qT "^(y ) X ^ T hv ) Xmnn,2 mnn,q(A. 2)For example, the angle between e and the first column of (A. 2) is givenXTcos m =^PT '^hv )mI, K'K^ ^,.yP = {p+l, p+2, . ..,p+q-1}hence.. 2=^"*p = Jp <V(\) \,i>^Jp 'Vf\' %,i' • ^ ''P- ^'orwhere'i' is of order 0(q_'^)P T>q^ = max ( | T (y ) | / | T (y ) | }, k ?f P^ ^ ' m k ' ' m p '11Therefore,<t>.^ is at most of order OCq.*^) |Proof of Theorem 2Taking the q vectors (A.l) each divided byT"^(vi), £ = p, p+1,.... p+q-1, as coordinate vectors W = [w , w ,^ , . . . , w . .] in E,, thep p+1 P+q-i J-2eigendirections of the projected operator A are the nxq matrix Y = WSwherein which.y"^ A"2 Y = Dj"^ (A. 3)D. = diagld » d , ..., d . ^ j,J P P+1 P+q-1Y^ Y = I,qand S is a qxq orthogonal matrix. ThusS*(wV^W)S = D."^,Jor(¥V^W)S = SD "^; (A.1+)t -2i.e., S is the eigenvector matrix of W A W, which can he written asW^A"^W = A"^ + K, (a. 5)where the elements of K are of order 0(t), i.e., 0[|t (y.)I/|T (y„)|]m 1 m X,2ji ^ P. Assume now A . = X. ^ = . . . = A, is an h-i+1 fold eigenvalue of1 1+1 hA; then as j-*^, h-i+1 independent eigensolutions of (A.i+) with d->X1exist. Everyone of these is described by q values s , s ,^, ..., s. ^,"p p+1 p+q-1and we assume that these solutions are normalized such that2 2 2s. +s., +...+S, =1. Then the s„ with £ ^ i, i+1, ..., h are of1 1+1 h Zh P'^g"!order 0(t). This means the angle between 2, ^o ^o ^-"^^ I ^o "^o ^^ alsoi ^ ^ ii=p ^ :hof order 0(t), while according to Theorem 1 the angle between J s wand the eigenspace of X., ^..-,> •••» X of A is of order 0(q. ). Thisestablishes the theorem. |12Appendix IIh =C C^8 Jwhere B and C are matrices of order 8,K.\ =a, 0.1k0.1 a, 0.10.1 0.10.1 ^Na,,C = diag(0.1+, . .. , O.U)and the centres of the Gerschgorin discs (a } are given "by(2, 3, 6, 9, 10, 11, 12, 13).^2 =16.778-U.889-1+.889-i;.889-1.556-0.222-li.88913.i+i+i^-1.556-1.5561.7783.111-1+.889-1.55613.i+i+i+-1.5561.7783.111-I1.889-1.556-1.55613.i+i+i+1.7783.111-1.5561.7781.7781.77810.1116Mh-0.2223.1113.1113.1116.1+1+1+8.778S-5 -h.1+ 61 -1+11+ 16 -1+-1+16-1+1,n = 61+13BIBLIOGRAPHIC DATASHEET1. Report No.UIUC-CAC-DN-73-913. Recipient's Accession No.4. Title and SubtitleON THE INTERMEDIATE EIGENVALUES OF SYMMETRICSPARSE MATRICES5. Report DateOctober 19736.7. Author(s)Ahmed Sameh, Jonathan Lermit, and Killion Noh8. Performine Organization Rept.Nc9. Performing Organization Name and AddressCenter for Advanced ComputationUniversity of Illinois at Urbana-ChampaignUrbana, Illinois 618OI10. Project/Task/Work Unit No.11. Contract/Grant No.DAHC04-72-C-000112. Sponsoring Organization Name and AddressU.S. Army Research OfficeDiike StationDurham, North Carolina13. Type of Report & PeriodCoveredResearch-interim14.15. Supplementary Notes16. AbstractsAn algorithm has been developed for finding the eigenvalues of asymmetric matrix A in a given interval [a, b] and the correspondingeigenvectors using a modification of the method of simultaneousiterations with the same favorable convergence properties. Thetechnique is most suitable for large sparse matrices and can be effectivelyimplemented on a parallel computer such as the ILLIAC IV.17. Key Words and Document Analysis. 17a. DescriptorsEigenvalues and EigenvectorsChebyshev PolynomialsSimultaneous Iterations17b. Identifiers/Open-Ended Terms17c. COSATI Field/Group18. Availability StatementORM NTIS-35 (REV. 3-72)No restriction on distribution.Available from National TechnicalInformation Service, Springfield, Va.221^119. Security Class (ThisReport)-^UNCI.ASSIFIFD20. Security Class (ThisPageUNCLASSIFIED21. No. of Pages22. PriceUSCOMM-DC 14952-P72

https://www.ideals.illinois.edu/bitstream/handle/2142/32233/onintermediateei91same.pdf?sequence=2

On the intermediate eigenvalues of symmetric sparse matrices, CAC No. 91

Abstract

Similar works

Full text

Available Versions

IDEALS @ Illinois

Illinois Digital Environment for Access to Learning and Scholarship Repository