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

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

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

. , n, is equal to the number of eigenvalues of A that are strictly greater than σ. Thisproperty allows the computation of the number of eigenvalues lying in a given interval.The determinants pk (σ) are especially easy to compute if A is tridiagonal, so A is usuallytransformed to this form before applying the Sturm sequence technique.4.4. GENERALIZED EIGENVALUE PROBLEMS4.4135Generalized Eigenvalue ProblemsMany eigenvalue problems occurring in practice have the form of a generalized eigenvalueproblemAx = λBx,where A and B are given n × n matrices.

In structural vibration problems, for example, Arepresents the stiffness matrix and B the mass matrix , and the eigenvalues and eigenvectorsdetermine the natural frequencies and modes of vibration of the structure (see ComputerProblem 4.12 for an example). A detailed study of the theory and algorithms for thisand other generalized eigenvalue problems is beyond the scope of this book, but the basicmethods available for their solution are briefly outlined next.If either of the matrices A or B is nonsingular, then the generalized eigenvalue problemcan be converted to a standard eigenvalue problem, either(B −1 A)x = λx or(A−1 B)x = (1/λ)x.Such a transformation is not generally recommended, however, since it may cause• Loss of accuracy due to rounding error in forming the product matrix, especially whenA or B is ill-conditioned• Loss of symmetry when A and B are symmetricIf A and B are symmetric, and one of them is positive definite, then symmetry can still beretained by using the Cholesky factorization.

For example, if B = LLT , then the generalized eigenvalue problem can be rewritten as the standard symmetric eigenvalue problem(L−1 AL−T )y = λy,and x can be recovered from the triangular linear system LT x = y. Transformation to astandard eigenvalue problem may still incur unnecessary rounding error, however, and itoffers no help if both A and B are singular.A numerically superior approach, which is applicable even when the matrices are singular or indefinite, is the QZ algorithm. Note that if A and B are both triangular, thenthe eigenvalues are given by λi = aii /bii for bii 6= 0. This circumstance is the motivationfor the QZ algorithm, which reduces A and B simultaneously to upper triangular form byorthogonal transformations.

First, B is reduced to upper triangular form by an orthogonaltransformation Q0 applied on the left, and the same orthogonal transformation is also applied to A. Then a sequence of orthogonal transformations Qk is applied to both matricesfrom the left to reduce A to upper Hessenberg form, and these alternate with orthogonaltransformations Zk applied on the right to restore B to upper triangular form. Finally, in aprocess analogous to QR iteration for the standard eigenvalue problem, additional orthogonal transformations are applied, alternating on the left and right, so that A converges toupper triangular form while maintaining the upper triangular form of B. The product ofall the transformations on the left is denoted by Q, and the product of those on the rightis denoted by Z, giving the algorithm its name.

The eigenvalues can now be determinedfrom the mutually triangular form, and the eigenvectors can be recovered via Q and Z.1364.54.5.1CHAPTER 4. EIGENVALUES AND SINGULAR VALUESSingular ValuesSingular Value DecompositionThe singular value decomposition (SVD) is an eigenvalue-like decomposition for rectangularmatrices. Let A be an m × n real matrix. Then the singular value decomposition has theformA = U ΣV T ,where U is an m × m orthogonal matrix, V is an n × n orthogonal matrix, and Σ is anm × n diagonal matrix, with0for i =6 jσij =.σi ≥ 0 for i = jThe diagonal entries σi are called the singular values of A and are usually ordered so thatσi ≥ σi+1 , i = 1, .

. . , n − 1. The columns ui of U and vi of V are the corresponding leftand right singular vectors.Example 4.16 Singular Value Decomposition. The singular value decomposition of1 2 3 4 5 6A= 7 8 910 11 12is given by U ΣV T =25.50.1410.825 −0.420 −0.351 0.344 00.4260.2980.782 0.5470.0280.664 −0.509   00.750 −0.371 −0.5420.079001.29000 0.5040  −0.76100.40800.574 0.644−0.057 0.646  .−0.816 0.408Thus, we have σ1 = 25.5, σ2 = 1.29, and σ3 = 0. A singular value of zero indicates that thematrix is rank-deficient; in general, the rank of a matrix is equal to the number of nonzerosingular values, which in this example is two.The singular values of A are the nonnegative square roots of the eigenvalues of AT A,and the columns of U and V are orthonormal eigenvectors of AAT and AT A, respectively.Algorithms for computing the SVD work directly with A, however, without forming AATor AT A, thereby avoiding any loss of information associated with forming these matrixproducts explicitly.The SVD is usually computed by a variant of QR iteration.

First, A is reduced tobidiagonal form by orthogonal transformations, then the remaining off-diagonal entries areannihilated iteratively. The SVD can also be computed by a variant of the Jacobi method,which can be useful on parallel computers or if the matrix has some special structure. Thetotal number of arithmetic operations required to compute the SVD of an m × n densematrix is proportional to mn2 + n3 , with the proportionality constants ranging from 2 to4.5.

SINGULAR VALUES13710 or more, depending on the particular algorithm used and the combination of singularvalues and right or left singular vectors desired. If the matrix is large and sparse, thenbidiagonalization is most effectively performed by a variant of the Lanczos algorithm, whichis especially suitable if only a few of the extreme singular values and corresponding singularvectors are needed.4.5.2Applications of SVDThe singular value decomposition A = U ΣV T has many important applications, amongwhich are the following:• Euclidean norm of a matrix.

The matrix norm subordinate to the Euclidean vector normis given by the largest singular value of the matrix,kAk2 = maxx6=okAxk2= σmax .kxk2• Condition number of a matrix. The condition number of a matrix A with respect to theEuclidean norm is given by the ratiocond(A) = σmax /σmin .This result agrees with the definition of cond(A) for a square matrix given in Section 2.3.3when using the Euclidean norm, and it also enables us to assign a condition number to arectangular matrix. Just as the condition number of a square matrix measures closenessto singularity, the condition number of a rectangular matrix measures closeness to rankdeficiency.• Rank of a matrix.

In theory, the rank of a matrix is equal to the number of nonzerosingular values it has. In practice, however, the rank may not be well-determined in thatsome singular values may be very small but nonzero. For many purposes it may be betterto regard any singular values falling below some threshold as negligible in determiningthe “numerical rank” of the matrix.

One way to interpret this is that the given matrix isvery near to (i.e., within the given threshold of) a matrix of the rank so determined.• Solving linear systems or linear least squares problems. The minimum Euclidean normsolution to Ax ≈ b is given byX uT bix=vi .σiσi 6=0The SVD is especially useful for ill-conditioned or rank-deficient problems, since “small”singular values can be dropped from the summation, thereby stabilizing the solution(making it much less sensitive to perturbations in the data).• Pseudoinverse of a matrix. Define the pseudoinverse of a scalar σ to be 1/σ if σ 6= 0,and zero otherwise.

Define the pseudoinverse of a (possibly rectangular) diagonal matrixby transposing the matrix and taking the scalar pseudoinverse of each entry. Then thepseudoinverse of a general real m × n matrix A, denoted by A+ , is given byA+ = V Σ+ U T .138CHAPTER 4. EIGENVALUES AND SINGULAR VALUESNote that the pseudoinverse always exists regardless of whether the matrix is square orof full rank. If A is square and nonsingular, then the pseudoinverse is the same as theusual matrix inverse, A−1 . In any case, the least squares solution to Ax ≈ b of minimumEuclidean norm is given by A+ b.• Orthonormal bases for range and null spaces.

The columns of V corresponding to zerosingular values form an orthonormal basis for the null space of A. The remaining columnsof V form an orthonormal basis for the orthogonal complement of the null space. Similarly, the columns of U corresponding to nonzero singular values form an orthonormalbasis for the range space of A, and the remaining columns of U form an orthonormalbasis for the orthogonal complement of the range space.• Approximating a matrix by one of lower rank. Another way to write the SVD isA = U ΣV T = σ1 E1 + σ2 E2 + · · · + σn En ,where Ei = ui viT . Each Ei is of rank 1 and can be stored using only m + n storagelocations.

Moreover, the product Ei x can be formed using only m + n multiplications.Thus, a useful condensed approximation to A can be obtained by omitting from theforegoing summation those terms corresponding to the smaller singular values, since theyhave relatively little effect on the sum. It can be shown that this approximation using thek largest singular values is the closest matrix of rank k to A in the Frobenius norm. (TheFrobenius norm of an m × n matrix is the Euclidean norm of the matrix considered as avector in Rmn .) Such an approximation is useful in image processing, data compression,cryptography, and numerous other applications.• Total least squares. In an ordinary linear least squares problem Ax ≈ b, we implicitlyassume that the entries of A are known exactly, whereas the entries of b are subjectto error.

In curve-fitting or regression problems where all of the variables are subjectto measurement error or other uncertainty, it may make more sense to minimize theorthogonal distances between the data points and the curve rather than the verticaldistances as in ordinary least squares. Such a total least squares solution can be computedusing the singular value decomposition [ A b ] = U ΣV T .

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

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