Nash - Compact Numerical Methods for Computers (523163), страница 26
Текст из файла (страница 26)
E is the diagonal matrixEij = ejδ i j(10.3)with δij the familiar Kronecker delta (see §6.2, p 60).By virtue of the orthogonality, we haveX T X = ln(10.4)where 1 n is the unit matrix of order n, but because X is composed of n orthogonaland non-null vectors, it is of full rank and invertible. Thus from equation (10.4)by left multiplication with X we getXXT X = X(10.5)-1so that right multiplication with X givesXXT XX-1 = XX T = XX -1 = 1 nT(10.6)-lshowing that X is the inverse X and that X is an orthogonal matrix.(iii) If the matrix A is not only symmetric, soAT = A(10.7)but also positive definite (see §7.1, p 71), then from the singular-value decompositionA = USVT(2.53)119Home Next120Compact numerical methods for computersthe eigenvalues of A are found as the diagonal elements Sii, i = 1, 2, .
. . , n, andthe matrices U and V both contain complete sets of eigenvectors. These sets (thecolumns of each matrix) are not necessarily identical since, if any two eigenvaluesare equal (also denoted as being degenerate or of multiplicity greater than one),any linear combination of their corresponding eigenvectors is also an eigenvectorof the same eigenvalue. If the original two eigenvectors are orthogonal, thenorthogonal linear combinations are easily generated by means of orthogonalmatrices such as the plane rotations of §3.3. This is an important point to keep inmind; recently a computer service bureau wasted much time and money trying tofind the ‘bug’ in a program taken from a Univac 1108 to an IBM 370 because theeigenvectors corresponding to degenerate eigenvalues were computed differentlyby the two computers.Property (iii) above will now be demonstrated.
First note that the eigenvalues ofa symmetric positive definite matrix are positive. For, from (7.9) and equation(10.1), we have0 < y TAy = yT XEXT y(10.8)thus implying that all the ej must be positive or else a vector w = XT y could bedevised such that wj = 1, wi = 0 for i j corresponding to the non-positiveeigenvalue, thus violating the condition (7.9) for definiteness.
Hereafter, E and Swill be ordered so thatSi > Si +l > 0(10.9)ei > ei + l ·(10.10)This enables S and E to be used interchangably once their equivalence has beendemonstrated.Now consider pre-multiplying equation (10.1) by A. Thus we obtainA2X = AAX = AXE = XEE = XE 2(10.11)while from symmetry (10.7) and the decomposition (2.53)A2 V = ATAV = VS 2 .(10.12)22Since (10.11) and (10.12) are both eigenvalue equations for A , S and E 2 areidentical to within ordering, and since all ei are positive, the orderings (10.9) and(10.10) imply(10.13)S = E.Now it is necessary to show that(10.14)AV = VS.From (10.1), letting Q = X T V, we obtainTAV = XEX V = XEQ = XSQ.However, from (10.11) and (10.12), we getQS2 = S2 Q.(10.15)(10.16)Real symmetric matrices121Explicit analysis of the elements of equation (10.16) shows that (a) if Sii Sjj, thenQij = 0, and (b) the commutationQS = SQ(10.17)is true even in the degenerate eigenvalue case; thus,AV = XSQ = XQS = XXTVS = VS.(10.18)The corresponding result for U is shown in the same fashion.10.2.
EXTENSION TO MATRICES WHICH ARE NOT POSITIVEDEFINITEThe above result shows that if A is a symmetric positive definite matrix, itseigensolutions can be found via the singular-value decomposition algorithms 1and 4. Moreover, only one of the matrices U or V is required, so that theeigenvectors overwrite the original matrix. Furthermore, the algorithm 1, forexample, generates the matrix B = US in performing an orthogonalisation of thecolumns of A, so that the eigenvalues are kept implicitly as the norms of thecolumns of the orthogonalised matrix and a separate vector to store the eigenvalues is not required.What of the case where the matrix A is not positive definite? This is hardly anyextra trouble, since the matrixA' = A – hl n(10.19)has the same eigenvectors as A (as proved by direct substitution in equation(10.1)) and has eigenvaluesE'ii = Eii – hfor i = 1, 2, .
. . , n(10.20)where Eii, i = 1, 2, . . . , n , are the eigenvalues of A. Thus it is always possible togenerate a positive definite matrix from A by adding an appropriate constant –hto each of its diagonal elements. Furthermore, it is very simple to compute anappropriate value for h from Gerschgorin’s theorem (Ralston 1965, p 467), whichstates that all the eigenvalues of a matrix A (a general square matrix, real orcomplex) are found in the domain which is the union of the circles having centresAii, i = 1, 2, .
. . , n , and respective radii(10.21)Because a symmetric matrix has real eigenvalues this impliesEnn > E = min(Aii – ri )(10.22)If E > 0, the matrix A is positive definite; otherwise a shift equal to E will make122Compact numerical methods for computersit so. Thus we can defineh=0h = –(|E| + ε½) = E – ε½for E > εfor E < ε(10.23a)(10.23b)to ensure a positive definite matrix A' results from the shift (10.19). The machineprecision ε is used simply to take care of those situations, such as a matrix with anull row (and column), where the lower bound E is in fact a small eigenvalue.Unfortunately, the accuracy of eigensolutions computed via this procedure issensitive to the shift.
For instance, the largest residual element R, that is, theelement of largest magnitude in the matrixAX – XE(10.24)and the largest inner product P, that is, the off-diagonal element of largestmagnitude in the matrix(10.25)XT X – 1 nfor the order-10 Ding Dong matrix (see appendix 1) are: for h = –3.57509,R = 5·36442E–6 and P = 1·24425E–6 while for h = –10·7238, R = 1·49012E–5and P = 2·16812E–6. These figures were computed on a Data General NOVA(23-bit binary arithmetic) using single-length arithmetic throughout as no extended precision was available.
The latter shift was obtained usingfor E > 0for E < 0.(10.26a)(10.26 b)In general, in a test employing all nine test matrices from appendix 1 of order 4and order 10, the shift defined by formulae (10.23) gave smaller residuals andinner products than the shift (10.26). The eigenvalues used in the above exampleswere computed via the Rayleigh quotient(10.27)rather than the singular value, that is, equation (10.20). In the tests mentionedabove, eigenvalues computed via the Rayleigh quotient gave smaller residualsthan those found merely by adding on the shift. This is hardly surprising if thenature of the problem is considered. Suppose that the true eigenvectors are φi ,i = 1, 2, .
. . , n. Let us add a component c w to φ j , where w is some normalisedcombination of the φ i , i j, and c measures the size of the component (error); thenormalised approximation to the eigenvector is thenxj = (l + c 2 ) - ½ ( φj + cw).(10.28)The norm of the deviation (xj – φ j ) is found, using the binomial expansion andignoring terms in c4 and higher powers relative to those in c2, to be approximatelyequal to c.
The Rayleigh quotient corresponding to the vector given by (10.28) is2 T2Qj = ( Ejj + c w Aw)/(1+ c )(10.29)sinceis zero by virtue of the orthogonality of the eigenvectors. Thedeviation of Qj from the eigenvalue isReal symmetric matrices123TABLE 10.1. Maximum absolute residual element R and maximum absolute inner product P betweennormalised eigenvectors for eigensolutions of order n = 10 real symmetric matrices. All programs in-22BASIC on a Data General NOVA. Machine precision = 2 .Algorithm 13 typeAlgorithm 14 typeRutishauserRutishauser with NashJacobiformulaeMatrixJacobiwhichordersJacobiusingsymmetrywithequation(10.27)withequation(10.20)HilbertRP7·26E-705·76E-68·64E-74·82E-61·13E-66·29E-61·10E-66·68E-62·32E-67·15E-62·32E-6Ding DongRP2·32E-602·86E-65·36E-78·86E-61·43E-61·08E-51·19E-65·36E-61·24E-61·54E-51·24E-6MolerRP1·74E-51·94E-73·62E-58·64E-76·34E-58·05E-71·01E-48·94E-73·91E-52·21E-69·46E-52·21E-6FrankRP2·29E-52·09E-75·53E-56·85E-78·96E-51·07E-61·25E-48·57E-75·72E-51·66E-69·72E-51·66E-6BorderedRP1·79E-65·34E-91·91E-65·96E-76·20E-69·98E-72·05E-51·40E-61·43E-65·54E-71·91E-65·54E-7DiagonalRP000000000000W+RP2·32E-61·79E-64·59E-61·26E-62·45E-51·88E-62·01E-51·91E-69·16E-61·75E-61·43E-51·75E-6W-RP1·94E-64·77E-78·58E-66·26E-71·63E-57·97E-72·86E-55·41E-71·35E-52·10E-62·00E-52·10E-6OnesRP4·65E-601·06E-63·65E-71·06E-59·92E-75·05E-51·04E-32·43E-59·89E-7l·l9E-59·89E-7(Ejj – Qj ) = c 2 (Ejj – w T Aw)/(1 + c2 ) .(10.30)Thus the error has been squared in the sense that the deviation of xj from φj is oforder c, while that of Qj from Ejj is of order c2.
Since c is less than unity, thisimplies that the Rayleigh quotient is in some way ‘closer’ to the eigenvalue thanthe vector is to an eigenvector.Unfortunately, to take advantage of the Rayleigh quotient (and residual calculation) it is necessary to keep a copy of the original matrix in the memory orperform some backing store manipulations. A comparison of results for algorithm13 using formulae (10.20) and (10.27) are given in table 10.1.Algorithm 13. Eigensolutions of a real symmetric matrix via the singular-valuedecompositionProcedure evsvd(n: integer; {order of matrix eigenproblem}var A,V : matrix; {matrix and eigenvectors}initev: boolean; {switch -- if TRUE eigenvectorsare initialized to a unit matrix of order n}124Compact numerical methods for computersAlgorithm 13.
Eigensolutions of a real symmetric matrix via the singular-valuedecomposition (cont.)W : wmatrix; {working array}var Z: rvector); {to return the eigenvalues}{alg13.pas ==eigensolutions of a real symmetric matrix via the singular valuedecomposition by shifting eigenvalues to form a positive definitematrix.This algorithm replaces Algorithm 13 in the first edition of CompactNumerical Methods.Copyright 1988 J.C.Nash}varcount, i, j, k, limit, skipped : integer;c, p, q, s, shift, t : real ; {rotation angle quantities}oki, okj, rotn : boolean;ch : char;beginwriteln(‘alg13.pas-- symmetric matrix eigensolutions via svd’);{Use Gerschgorin disks to approximate the shift.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
