Nash - Compact Numerical Methods for Computers (523163), страница 23
Текст из файла (страница 23)
Presuming e1 is the largesteigenvalue in magnitude, xi+1 can also be written(9.6)But since| ej/e1 | < 1(9.7)unless j = 1 (the case of degenerate eigenvalues is treated below), the coefficientsof φj, j 1, eventually become very small. The ultimate rate of convergence isgiven byr = | e 2 /e1 |(9.8)where e 2 is the eigenvalue having second largest magnitude. By working with thematrixA ' = A – kl(9.9)this rate of convergence can be improved if some estimates of e 1 and e2 areknown. Even if such information is not at hand, ad hoc shifts may be observed toimprove convergence and can be used to advantage.
Furthermore, shifts permit(i) the selection of the most positive eigenvalue or the most negative eigenvalueand, in particular,(ii) evasion of difficulties when these two eigenvalues are equal in magnitude.Degenerate eigenvalues present no difficulty to the power method except that itnow converges to a vector in the subspace spanned by all eigenvectors corresponding to e1 . Specific symmetry or other requirements on the eigenvector mustbe imposed separately.In the above discussion the possibility that a 1 = 0 in the expansion of x 1 hasbeen conveniently ignored, that is, some component of x 1 in the direction of φ l i sassumed to exist.
The usual advice to users is, ‘Don’t worry, rounding errors will104Compact numerical methods for computerseventually introduce a component in the right direction’. However, if the matrixA has elements which can be represented exactly within the machine, that is, if Acan be scaled so that all elements are integers small enough to fit in one machineword, it is quite likely that rounding errors in the ‘right direction’ will not occur.Certainly such matrices arise often enough to warrant caution in choosing astarting vector.
Acton (1970) and Ralston (1965) discuss the power method inmore detail.The power method is a simple, yet highly effective, tool for finding the extremeeigensolutions of a matrix. However, by applying it with the inverse of the shiftedmatrix A' (9.9) an algorithm is obtained which permits all distinct eigensolutionsto be determined.
The iteration does not, of course, use an explicit inverse, butsolves the linear equationsA 'y i = xi(9.10a)xi+l = yi/||yi||.(9.10b)then normalises the solution byNote that the solution of a set of simultaneous linear equations must be found ateach iteration.While the power method is only applicable to the matrix eigenproblem (2.62),inverse iteration is useful for solving the generalised eigenproblem (2.63) usingA' = A – kB(9.11)in place of (9.9). The iteration scheme is nowA' yi = Bxi(9.12a)xi+1 = yi/||yi||.(9.12b)Once again, the purpose of the normalisation of y in (9.1b), (9.10b) and (9.12b) issimply to prevent overflow in subsequent calculations (9.1a), (9.10a) or (9.12a).The end use of the eigenvector must determine the way in which it is standardised. In particular, for the generalised eigenproblem (2.63), it is likely that xshould be normalised so thatx T Bx = 1.(9.13)Such a calculation is quite tedious at each iteration and should not be performeduntil convergence has been obtained, since a much simpler norm will suffice, forinstance the infinity norm(9.14)where yj is the jth element of y.
On convergence of the algorithm, the eigenvalueis(9.14)e = k + xj/y j(where the absolute value is not used).Inverse iteration works by the following mechanism. Once again expand x 1 asThe algebraic eigenvalue problem105in (9.2); then(9.16)or(9.17)Therefore(9.18)and the eigenvector(s) corresponding to the eigenvalue closest to k very quicklydominate(s) the expansion. Indeed, if k is an eigenvalue, A' is singular, and aftersolution of the linear equations (9.12a) (this can be forced to override thesingularity) the coefficient of the eigenvector φ corresponding to k should be ofthe order of 1/eps, where eps is the machine precision.
Peters and Wilkinson(1971, pp 418-20) show this ‘full growth’ to be the only reliable criterion forconvergence in the case of non-symmetric matrices. The process then convergesin one step and obtaining full growth implies the component of the eigenvector inthe expansion (9.2) of x 1 is not too small.
Wilkinson proposes choosing differentvectors x 1 until one gives full growth. The program code to accomplish this isquite involved, and for symmetric matrices repetition of the iterative step issimpler and, because of the nature of the symmetric matrix eigenproblem, canalso be shown to be safe.
The caution concerning the choice of starting vector formatrices which are exactly representable should still be heeded, however. In thecase where k is not an eigenvalue, inverse iteration cannot be expected toconverge in one step. The algorithm given below therefore iterates until thevector x has converged.The form of equation (9.12a) is amenable to transformation to simplify theiteration. That is, pre-multiplication by a (non-singular) matrix Q givesQAyi = QBxi.(9.19)Note that neither xi nor yi are affected by this transformation. If Q is taken to bethe matrix which accomplishes one of the decompositions of §2.5 then it isstraightforward to carry out the iteration.
The Gauss elimination, algorithm 5, orthe Givens’ reduction, algorithm 3, can be used to effect the transformation forthis purpose. The matrix Q never appears explicitly. In practice the two matricesA and B will be stored in a single working array W. Each iteration will correspondto a back-substitution similar to algorithm 6. One detail which deserves attentionis the handling of zero diagonal elements in QA', since such zeros imply a divisionby zero during the back-substitution which solves (9.19). The author has foundthat replacement of zero (or very small) elements by some small number, say themachine precision multiplied by the norm of A', permits the process to continue106Compact numerical methods for computersquite satisfactorily. Indeed I have performed a number of successful computationson a non-symmetric generalised eigenvalue problem wherein both A and B weresingular by means of this artifice! However, it must be admitted that thisparticular problem arose in a context in which a great deal was known about theproperties of the solutions.Algorithm 10.
Inverse iteration via Gauss eliminationThe algorithm requires a working array W, n by 2 * n and two vectors x and y of order n.procedure gii(nRow : integer; {order of problem}var A : rmatrix; {matrices of interest packedinto array as ( A | B) }var Y : rvector; {the eigenvector}var shift : real; {the shift -- on input, the valuenearest which an eigenvalue is wanted. On output, theeigenvalue approximation found.}var itcount: integer); {On input a limit to the numberof iterations allowed to find the eigensolution.
Onoutput, the number of iterations used. The returnedvalue is negative if the limit is exceeded.}{alg10.pas == Inverse iteration to find matrix an eigensolution ofAY=ev*B*Yfor real matrices A and B of order n. The solution found corresponds toone of the eigensolutions having eigenvalue, ev, closest to the valueshift. Y on input contains a starting vector approximation.Copyright 1988 J.C.Nash}vari, itlimit, j, k, m, msame, nRHS :integer;ev, s, t, to1 : real; {eigenvalue approximation}X : rvector;beginitlimit:=itcount; {to get the iteration limit from the call}nRHS:=nRow; {the number of right hand sides is also n since we willstore matrix B in array A in the right hand nRow columns}tol:=Calceps;s:=0.0; {to initialize a matrix norm}for i:=1 to nRow dobeginX[i]:=Y[i]; {copy initial vector to X for starting iteration}Y[i]:=0.0; {to initialize iteration vector Y}for j:=1 to nRow dobeginA[i,j]:=A[i,j]-shift*A[i,j+nRow];s:=s+abs(A[i,j]);end;end;tol:=tol*s; {to set a reasonable tolerance for zero pivots}gelim(nRow, nRHS, A, tol); {Gauss elimination STEPS 2-10}itcount:=0,msame :=0; {msame counts the number of eigenvector elements whichare unchanged since the last iteration}The algebraic eigenvalue problemAlgorithm 10.
Inverse iteration via Gauss elimination (cont.)while (msame<nRow) and (itcount<itlimit) dobegin {STEP 11 -- perform the back-substitution first}itcount:=itcount+1; {to count the iterations}m:=nRow; s:=X[nRow];X[nRow]:=Y[nRow]; {save last trial vector -- zeros on iteration 1}if abs(A[nRow,nRow])<tol then Y[nRow]:=s/tolelse Y[nRow]:=s/A[nRow,nRow];t:=abs(Y[nRow]);{ to set the first trial value for vector Y}for i:=(nRow-1) downto 1 do {STEP 12}begin {back-substition}s:=X[i]; X[i]:=Y[i];for j:=(i+1) to nRow dobegins:=s-A[i,j]*Y[j];end;if abs(A[i,i])<tol then Y[i]:=s/tol else Y[i]:=s/A[i,i];if abs(Y[i])>t thenbeginm:=i; t:=abs(Y[i]);end; {to update new norm and its position}end; {loop on i}ev:=shift+X[m]/Y[m]; {current eigenvalue approximation -- STEP 13}writeln(‘Iteration’ ,itcount,’ approx. ev=’,ev);{Normalisation and convergence tests -- STEP 14}t:=Y[m]; msame:=0;for i:=1 to nRow dobeginY[i]:=Y[i]/t;if reltest+Y[i] = reltest+X[i] then msame:=msame+1;{This test is designed to be machine independent.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
