Heath - Scientific Computing (523150), страница 35

Просмтор этого файла доступен только зарегистрированным пользователям. Но у нас супер быстрая регистрация: достаточно только электронной почты!

Текст из файла (страница 35)

.. ......  000···1 −a0 −a1 −a2 · · · −an−1using the methods discussed in this chapter.Although it is not useful numerically, the characteristic polynomial does permit us tomake an important theoretical observation about computing eigenvalues. Abel proved thatthe roots of a polynomial of degree greater than four cannot always be expressed by a closedform formula in the coefficients using ordinary arithmetic operations and root extractions.Thus, in general, computing the eigenvalues of matrices of order greater than four requiresa (theoretically infinite) iterative process.Example 4.5 Characteristic Polynomial.

To illustrate some of the numerical difficulties associated with the characteristic polynomial, consider the matrix1 A=, 1where is a positive number slightly smaller than the square root of machine precision ina given floating-point system. The exact eigenvalues of A are 1 + and 1 − . Computingthe characteristic polynomial of A in floating-point arithmetic, we getdet(A − λI) = λ2 − 2λ + (1 − 2 ) = λ2 − 2λ + 1,which has 1 as a double root. Thus, we cannot resolve the two eigenvalues by this methodeven though they are quite distinct in the working precision.

We would need up to twicethe precision in the coefficients of the characteristic polynomial to compute the eigenvaluesto the same precision as that of the input matrix.1224.2.2CHAPTER 4. EIGENVALUES AND SINGULAR VALUESJacobi Method for Symmetric MatricesOne of the oldest methods for computing eigenvalues of symmetric matrices is due to Jacobi.Starting with a symmetric matrix A0 = A, each iteration has the formAk+1 = JkT Ak Jk ,where Jk is a plane rotation chosen to annihilate a symmetric pair of entries in the matrixAk (so that the symmetry of the original matrix is preserved). Recall from Section 3.4.5that a plane rotation is an orthogonal matrix that differs from the identity matrix I in onlyfour entries, and this 2 × 2 submatrix has the formc s,−s cwith c and s the cosine and sine of the angle of rotation, respectively, so that c2 + s2 = 1.The choice of c and s is slightly more complicated in this context than in the Givensmethod for QR factorization because we are annihilating a symmetric pair of matrix entriesby a similarity transformation, as opposed to annihilating a single entry by a one-sidedtransformation.As before, it suffices to consider only the 2 × 2 case,c −sa bc sTJ AJ =scb d−s cc2 a − 2csb + s2 dc2 b + cs(a − d) − s2 b,=c2 b + cs(a − d) − s2 bc2 d + 2csb + s2 awhere b 6= 0 (else there is nothing to do).

The transformed matrix will be diagonal ifc2 b + cs(a − d) − s2 b = 0.Dividing both sides of this equation by c2 b, we obtain1+s (a − d) s2− 2 = 0.cbcMaking the substitution t = s/c, we obtain a quadratic equation1+t(a − d)− t2 = 0b√for t, the tangent of the angle of rotation, from which we can recover c = 1/ 1 + t2 ands = c · t. It is advantageous numerically to use the root of smaller magnitude of the equationfor t.Example 4.6 Plane Rotation. To illustrate the use of a plane rotation to annihilate asymmetric pair of off-diagonal entries, we consider the 2 × 2 matrix1 2A=.2 14.2. METHODS FOR COMPUTING ALL EIGENVALUES123The quadratic equation for the tangent reduces to t2 = 1 in this case, so we have t = ±1.Since the√ two roots are√of the same magnitude, we arbitrarily choose t = −1, which yieldsc = 1/ 2 and s = −1/ 2.

Using the resulting plane rotation J , we then have√√ √ √1/√2 1/√21 21/√2 −1/√230TJ AJ ==.−1/ 2 1/ 22 11/ 21/ 20 −1In the Jacobi method, plane rotations determined in this manner are repeatedly applied from both sides in systematic sweeps through the matrix until the off-diagonal massof the matrix is reduced to within some tolerance of zero. The resulting approximatelydiagonal matrix is orthogonally similar to the original matrix; hence, we have the approximate eigenvalues on the diagonal, and the product of all of the plane rotations gives theeigenvectors.Although the Jacobi method is reliable, simple to program, and capable of very highaccuracy, it converges rather slowly. It is also difficult to generalize beyond symmetric (orHermitian) matrices.

Except for very small problems, the Jacobi method usually requiresfive to ten times more work than more modern methods. Recently, however, the Jacobimethod has regained popularity because it is relatively easy to implement on parallel computers.The main source of inefficiency in the Jacobi method is that entries that have been annihilated by a previous iteration can subsequently become nonzero again, thereby requiringrepeated annihilation. The main computational advantage of more modern methods is thatthey are carefully designed to preserve zero entries once they have been introduced into thematrix.Example 4.7 Jacobi Method.

Let1 0A0 = 0 22 121.1We will repeatedly sweep through the matrix by rows and columns, annihilating successivematrix entries. We first annihilate the symmetrically placed entries (1, 3) and (3, 1) usingthe plane rotation0.707 0 −0.70730.7070J0 =  010  to obtain A1 = J0T A0 J0 =  0.70720.707  .0.707 00.70700.707 −1We next annihilate the symmetrically placed entries (1, 2) and (2, 1) using the plane rotation0.888 −0.460 03.36600.325J1 =  0.4600.888 0  to obtain A2 = J1T A1 J1 =  01.634 0.628  .0010.325 0.628−1We next annihilate the symmetrically placed entries (2, 3) and (3, 2) using the plane rotation1003.366 0.0720.317J2 =  0 0.975 −0.221  to obtain A3 = J2T A2 J2 =  0.072 1.7760 .0 0.2210.9750.3170−1.142124CHAPTER 4.

EIGENVALUES AND SINGULAR VALUESBeginning a new sweep, we again annihilate the symmetrically placed1) using the plane rotation0.998 0 −0.0703.388J3 =  010  to obtain A4 = J3T A3 J3 =  0.0720.070 00.9980entries (1, 3) and (3,0.0721.7760.00500.005  .−1.164This process continues until the off-diagonal entries are reduced to as small a magnitude asdesired. The result is an approximately diagonal matrix that is orthogonally similar to theoriginal matrix, with the orthogonal similarity transformation given by the product of theplane rotations.4.2.3QR IterationQR iteration for computing eigenvalues and eigenvectors makes repeated use of the QRfactorization to produce a unitary similarity transformation of the matrix to diagonal ortriangular form.

We initially take A0 = A, then at iteration k we compute the QR factorizationAk = Qk Rk ,and then form the reverse productAk+1 = Rk Qk .SinceRk Qk = QHk Ak Qk ,we see that the successive matrices Ak are unitarily similar to each other. Moreover, it canbe shown that if the moduli of the eigenvalues are all distinct, then the Ak converge totriangular form for a general initial matrix, or diagonal form for a symmetric initial matrix.This condition is not a serious restriction in practice because eigenvalues are deflated outas they are determined, which in this context simply means that we deal with successivelysmaller submatrices as each eigenvalue is determined. For example, after the last rowof the matrix has converged to within some tolerance of zero (except, of course, for thediagonal entry, which is then an approximate eigenvalue), it need be processed no furtherand attention can be restricted to the leading submatrix of dimension n−1.

Such reductionscontinue successively until all the eigenvalues have been found.Example 4.8 QR Iteration. Let7 2A0 =.2 4We first compute the QR factorization A0 = Q0 R0 , obtaining0.962 −0.2757.28Q0 =and R0 =0.2750.96203.02.3.304.2. METHODS FOR COMPUTING ALL EIGENVALUES125We next form the reverse product7.83A1 = R0 Q0 =0.9060.906.3.17We see that the off-diagonal entries are now smaller, so that the matrix is closer to beingdiagonal, and the diagonal entries are now closer to the eigenvalues, which are 8 and 3 forthis problem.

Repetition of this process would continue until the matrix is within toleranceof being diagonal, and the diagonal entries would then closely approximate the eigenvalues.The product of the orthogonal matrices Qk would yield the corresponding eigenvectors.The convergence rate of QR iteration can be accelerated by incorporating shifts of thefollowing form:Ak − σk I = Qk Rk ,Ak+1 = Rk Qk + σk I,where σk is a rough approximation to an eigenvalue.

This is called a shift because the entirespectrum of the matrix is displaced temporarily by the amount σk and then subsequentlyrestored. One choice for the shift is simply the lower right corner entry of the matrix. Abetter shift can be determined by computing the eigenvalues of the 2 × 2 submatrix in thelower right corner of the matrix. In either case, such a shift will become increasingly better(i.e., closer to an eigenvalue) as the matrix converges to diagonal or triangular form.Example 4.9 QR Iteration with Shifts. To illustrate the QR algorithm with shifts, werepeat the previous example using a shift of σ0 = 4, which is the lower right corner entryof the matrix.

Thus, we first compute the QR factorization A0 − σ0 I = Q0 R0 so that wehave0.8320.5553.61 1.66Q0 =and R0 =.0.555 −0.83201.11We next form the reverse product and add back the shift to obtain7.92 0.615A1 = R0 Q0 + σ0 I =.0.615 3.08Compared with the unshifted algorithm, the off-diagonal entries are smaller after one iteration, and the diagonal entries are closer approximations to the eigenvalues.

For the nextiteration, we would use the new value of the lower right corner entry as the shift.4.2.4Preliminary ReductionIn the simple form just given, each iteration of the QR method requires O(n3 ) work. Thework per iteration can be reduced if the matrix is initially transformed into a simpler form.In particular, it is advantageous if the matrix is as close as possible to triangular (or diagonalfor a symmetric matrix) before the QR iterations begin. A Hessenberg matrix is triangularexcept for one additional nonzero diagonal immediately adjacent to the main diagonal.Note that a symmetric Hessenberg matrix is tridiagonal.

Any matrix can be reduced to126CHAPTER 4. EIGENVALUES AND SINGULAR VALUESHessenberg form in a finite number of steps by an orthogonal similarity transformation, forexample, by using Householder transformations. Moreover, upper Hessenberg or tridiagonalform can then be preserved during the subsequent QR iterations. The advantages of thisinitial reduction to upper Hessenberg or tridiagonal form are• The work per QR iteration is reduced to at most O(n2 ).• The convergence rate of the QR iterations is enhanced.• If there are any zero entries on the first subdiagonal, then the problem can be brokeninto two or more smaller subproblems.Thus, the QR method is usually implemented as a two-stage process:Symmetric−→General−→tridiagonalorHessenberg−→diagonal−→triangularThe preliminary reduction requires a definite number of steps, whereas the subsequentiterative stage continues until convergence.

In practice, however, only a modest number ofiterations is usually required, so the O(n3 ) cost of the preliminary reduction is a significantfraction of the total. The total cost is strongly affected by whether the eigenvectors areneeded because their inclusion determines whether the orthogonal transformations must beaccumulated. For the symmetric case, the overall cost is roughly 34 n3 arithmetic operations(counting both additions and multiplications) if only the eigenvalues are needed, and about9n3 operations if the eigenvectors are also desired.

Характеристики

Список файлов книги