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

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

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

Experiment with variousvalues for the scaling parameter α. How do theaccuracy and execution time of this methodcompare with those of the normal equationsand QR factorization methods?3.8 The covariance matrix for the m × nleast squares problem Ax ≈ b is given byσ 2 (AT A)−1 , where σ 2 = kb − Axk22 /(m − n)at the least squares solution x. The entriesof this matrix contain important informationabout the goodness of the fit and any crosscorrelations among the fitted parameters. Thecovariance matrix is an exception to the general rule that inverses of matrices should neverbe computed explicitly. If an orthogonalization method is used to solve the least squaresproblem, then the normal equations matrixAT A is never formed, so we need an alternative method for computing the covariancematrix.(a) Show that (AT A)−1 = (RT R)−1 , whereR is the upper triangular factor obtained byQR factorization of A.(b) Based on this fact, implement a routine forcomputing the covariance matrix using onlythe already computed R.

(For purposes of thisexercise, you may ignore the scalar factor σ 2 .)Test your routine on a few example matricesto confirm that it gives the same result as computing (AT A)−1 .3.9 Most library routines for computing theQR factorization of an m × n matrix A return the matrix R in the upper triangle of thestorage for A and the Householder vectors inthe lower triangle of A, with an extra vectorto accommodate the overlap on the diagonal.Write a routine that takes this output arrayand auxiliary vector and forms the orthogonalmatrix Q explicitly by multiplying the corresponding sequence of Householder transformations times an m × m matrix that is initializedto the identity matrix I. Of course, the latter will require a separate array.

Test yourprogram on several randomly chosen matricesand confirm that your computed Q is indeedorthogonal and that the product114CHAPTER 3. LINEAR LEAST SQUARESQROrecovers A.3.10 (a) Implement both the classical andmodified Gram-Schmidt procedures and useeach to generate an orthogonal matrix Qwhose columns form an orthogonal basis forthe column space of the Hilbert matrix H,with entries hij = 1/(i + j − 1), for n =2, . .

. , 12 (see Computer Problem 2.6). As ameasure of the quality of the results (specifically, the potential loss of orthogonality), plotthe quantity − log10 (kI − QT Qk), which canbe interpreted as “digits of accuracy,” for eachmethod as a function of n.

In addition, tryapplying the classical procedure twice (i.e., apply your classical Gram-Schmidt routine to itsown output Q to obtain a new Q), and againplot the resulting departure from orthogonality. How do the three methods compare inspeed, storage, and accuracy?(b) Repeat the previous experiment, but thistime use the Householder method, that is, usethe explicitly computed orthogonal matrix Qresulting from Householder QR factorizationof the Hilbert matrix. Note that if the routine you use for Householder QR factorizationdoes not form Q explicitly, then you can obtain Q by multiplying the sequence of Householder transformations times a matrix that isinitialized to the identity matrix I (see previous exercise).

Again, plot the departure fromorthogonality for this method and compare itwith that of the previous methods.(c) Yet another way to compute an orthogonalbasis is to use the normal equations. If we formthe normal equations matrix and compute itsCholesky factorization AT A = LLT , then wehaveI= L−1 (AT A)L−T=(AL−T )T (AL−T ),which means that Q = AL−T is orthogonal,and its column space is obviously the same asthat of A. Repeat the previous experiment using Hilbert matrices again, this time using theQ obtained in this way from the normal equations (the required triangular solution may bea little tricky, depending on the software youuse).

Again, plot the resulting departure fromorthogonality and compare it with that of theprevious methods.(d ) Can you explain the relative quality of theresults you obtained for the various methodsused in these experiments?3.11 What is the exact solution to the linearleast squares problem 11 1 1  x1 0 0  0 x2 ≈  0 00x30 0 0as a function of ?Solve this least squares problem using each ofthe following methods. For each method, experiment with the value of the parameter tosee how small you can take it and still obtainan accurate solution.

Pay particular attention√to values around ≈ mach and ≈ mach .(a) Normal equations method(b) Augmented system method(c) Householder QR method(d ) Givens QR method(e) Classical Gram-Schmidt orthogonalization(f ) Modified Gram-Schmidt orthogonalization(g) Classical Gram-Schmidt orthogonalizationwith iterative refinement (i.e., CGS appliedtwice)Chapter 4Eigenvalues and Singular Values4.1Eigenvalues and EigenvectorsThe standard algebraic eigenvalue problem is as follows: Given an n × n matrix A, find ascalar λ and a nonzero vector x such thatAx = λx.Such a scalar λ is called an eigenvalue, and x is a corresponding eigenvector .

In addition tothe “right” eigenvector defined above, we could also define a “left” eigenvector y such thaty T A = λy T , but since a left eigenvector of A is a right eigenvector of AT , we will consideronly right eigenvectors. The set of all the eigenvalues of a matrix A, denoted by λ(A), iscalled the spectrum of A. The maximum modulus of the eigenvalues, max{|λ|: λ ∈ λ(A)},is called the spectral radius of A, denoted by ρ(A).An eigenvector of a matrix determines a direction in which the effect of the matrix isparticularly simple: The matrix expands or shrinks any vector lying in that direction bya scalar multiple, and the expansion or contraction factor is given by the correspondingeigenvalue λ. Thus, eigenvalues and eigenvectors provide a means of understanding thecomplicated behavior of a general linear transformation by decomposing it into simpleractions.Eigenvalue problems occur in many areas of science and engineering.

For example, thenatural modes and frequencies of vibration of a structure are determined by the eigenvectors and eigenvalues of an appropriate matrix. The stability of the structure is determinedby the locations of the eigenvalues, and thus their computation is of critical interest. Wewill also see later in this book that eigenvalues can be very useful in analyzing numericalmethods, such as the convergence analysis of iterative methods for solving systems of algebraic equations, and the stability analysis of methods for solving systems of differentialequations.Although most of our examples will involve only real matrices, both the theory andcomputational procedures we will discuss in this chapter are generally applicable to complexmatrices.

Notationally, the only difference in dealing with complex matrices is that the115116CHAPTER 4. EIGENVALUES AND SINGULAR VALUESconjugate transpose, denoted by AH , is used instead of the usual matrix transpose, AT(recall the definitions of transpose and conjugate transpose from Section 2.5).Example 4.1 Eigenvalues and Eigenvectors. 01 01.1. A =: λ = 1, x =and λ = 2, x =10 20 11 11.2.

A =: λ = 1, x =and λ = 2, x =010 2 3 −1113. A =: λ = 2, x =and λ = 4, x =.−131−1 1.5 0.51−14. A =: λ = 2, x =and λ = 1, x =.0.5 1.511 √0 11i5. A =: λ = i, x =and λ = −i, x =, where i = −1.−1 0i1Note that for examples 1 and 2 the eigenvalues are the diagonal entries of A, and for example1 the eigenvectors are the columns of the identity matrix I. The matrices in examples 3and 4 are symmetric, and the eigenvalues are real. Example 5 shows, however, that anonsymmetric real matrix need not have real eigenvalues.4.1.1NonuniquenessNeither eigenvalues nor eigenvectors are necessarily unique, in the following senses:• The eigenvalues of a matrix are not necessarily all distinct.

That is, more than onedirection may have the same expansion or contraction factor. In this case, we say thatthe matrix has a multiple eigenvalue. For example, 1 is an eigenvalue of multiplicity nfor the n × n identity matrix I.• Eigenvectors can obviously be scaled arbitrarily: if Ax = λx, then A(γx) = λ(γx) forany scalar γ, so that γx is also an eigenvector corresponding to λ. For example, 1 11γIf A =, then γx = γ=0 21γis an eigenvector corresponding to the eigenvalue λ = 2 for any nonzero scalar γ.

Consequently, eigenvectors are usually normalized by requiring some norm of the vector to be1.4.1.2Characteristic PolynomialThe equation Ax = λx is equivalent to(A − λI)x = o.4.1. EIGENVALUES AND EIGENVECTORS117This homogeneous equation has a nonzero solution x if and only if its matrix is singular.Thus, the eigenvalues of A are the values λ such thatdet(A − λI) = 0.Now det(A − λI) is a polynomial of degree n in λ, called the characteristic polynomial ofA, and its roots are the eigenvalues of A.Example 4.2 Characteristic Polynomial. As an example, consider the characteristicpolynomial of one of the matrices in Example 4.1:3 −11 03−λ−1det−λ= det0 1−13−λ−13= (3 − λ)(3 − λ) − (−1)(−1) = λ2 − 6λ + 8 = 0,so that the eigenvalues are given by√6 ± 36 − 326±2λ===222and 4.Because the eigenvalues of a matrix are the roots of its characteristic polynomial, wecan conclude from the Fundamental Theorem of Algebra that an n × n matrix A alwayshas n eigenvalues, but they need be neither distinct nor real.

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

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