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

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

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

Although the modified Gram-Schmidt procedure has advantages in somecircumstances, for solving least squares problems it is somewhat inferior to the Householdermethod in storage, work, and accuracy.Example 3.8 Gram-Schmidt QR Factorization. We illustrate modified Gram-Schmidtorthogonalization by again solving the quadratic polynomial data-fitting problem in Example 3.2, with1 −1.0 1.01.0 1 −0.5 0.25  0.5 A = 10.0 0.0  , b =  0.0  .1 0.5 0.5 0.25 11.0 1.02.0Normalizing the first column of A, we computer1,1 = ka1 k2 = 2.236,q1 = a1 /r1,10.447 0.447 = 0.447  . 0.447 0.447Orthogonalizing the first column against the subsequent columns, we getr1,2 = q1T a2 = 0,r1,3 = q1T a3 = 1.118,so that the matrix is transformed to become0.477 −1.0 0.447 −0.5 0.4470.0 0.4470.50.4471.00.50−0.25 −0.50 .−0.25 0.50Normalizing the second column, we compute−0.632 −0.316 .0= 0.316 0.632r2,2 = ka2 k2 = 1.581,q2 = a2 /r2,2Orthogonalizing the second column against the third column, we getr2,3 = q2T a3 = 0,3.4.

ORTHOGONALIZATION METHODS101so that the third column is unaffected. Finally, we normalize the third column0.535 −0.267 r3,3 = ka3 k2 = 0.935, q3 = a3 /r3,3 =  −0.535  . −0.267 0.535We have thus obtained the QR factorization0.447 −0.6320.535  0.447 −0.316 −0.267  2.236 0A=0−0.535  0.447 0.4470.316 −0.267 00.4470.6320.53501.58101.1180  = QR.0.935The transformed right-hand side is obtained from0.4470.447Q b = −0.632 −0.3160.535 −0.267T0.4470−0.5350.4470.316−0.267 1.00.447 1.789 0.5  0.632   0.0  = 0.632 .0.5350.51.3362.0We can now solve the upper triangular system Rx = QT b by back-substitution to obtain0.086x =  0.400  .1.4293.4.7Rank DeficiencySo far we have assumed that A is of full rank, rank(A) = n. If this is not the case, i.e.,if A has linearly dependent columns, then the QR factorization still exists, but the uppertriangular factor R is singular (as is AT A).

Thus, many vectors x give the same minimumresidual norm, and the least squares solution is not unique. This situation usually arisesfrom a poorly designed experiment, insufficient data, or an inadequate or redundant model.Thus, the problem should probably be reformulated or rethought.If one insists on forging ahead as is, however, then a common practice is to selectthe minimum residual solution x having the smallest norm. This may be computed byQR factorization with column pivoting, which we consider next, or by the singular valuedecomposition (SVD), which we will study in Section 4.5. Note that such a procedure fordealing with rank deficiency also enables us to handle underdetermined problems, wherem < n, since the columns of A are necessarily linearly dependent in that case.In practice, the rank of a matrix is often not clear-cut. Thus, a relative tolerance isused to detect near rank deficiency of least squares problems, just as in detecting near102CHAPTER 3.

LINEAR LEAST SQUARESsingularity of square linear systems. If a least squares problem is nearly rank-deficient,then the solution will be sensitive to perturbations in the input data. We will be able toexamine these issues more precisely when we introduce the singular value decomposition ofa matrix in Section 4.5. Within the context of QR factorization, the most robust method fordetecting and dealing with possible rank deficiency is column pivoting, which we considernext.Example 3.9 Near Rank Deficiency.

Consider the 3 × 2 matrix0.641 0.242A =  0.321 0.121  .0.962 0.363If we compute the QR factorization of A, we find that1.1997 0.4527R=.00.0002Thus, R is extremely close to being singular (indeed, it is exactly singular to the threedigit accuracy with which the problem was stated), and if we use R to solve a least squaresproblem, the result will be correspondingly sensitive to perturbations in the right-hand side.For practical purposes, the rank of A is only one rather than two, since its columns arenearly linearly dependent.3.4.8Column PivotingThe columns of a matrix A can be viewed as an unordered set of vectors from which we wishto select a maximal linearly independent subset.

Rather than processing the columns in thenatural order in computing the QR factorization, we instead select for reduction at eachstage the column of the remaining unreduced submatrix having maximum Euclidean norm.This column is interchanged (explicitly or implicitly) with the next column in the naturalorder and then is zeroed below the diagonal in the usual manner. The transformationrequired to do this must then be applied to the remaining unreduced columns, and theprocess is repeated. The process just described is called column pivoting.If rank(A) = k < n, then after k steps of this procedure, the norms of the remainingunreduced columns will be zero (or “negligible” in finite-precision arithmetic) below row k.Thus, we have produced an orthogonal factorization of the formR STQ AP =,O Owhere R is k × k, upper triangular, and nonsingular, and P is a permutation matrix thatperforms the column interchanges.

At this point, a basic solution (i.e., a solution havingat most k nonzero components) to the least squares problem Ax ≈ b can be computedby solving the triangular system Ry = c, where c is a vector composed of the first kcomponents of QT b, and then taking yx=P.o3.5. COMPARISON OF METHODS103In the context of data fitting, this procedure amounts to ignoring components of the modelthat are redundant or not well-determined.

If a minimum-norm solution is desired, however, it can be computed at the expense of some additional processing (from the right) toannihilate S as well.In practice, the rank of A is usually unknown, so the column pivoting process is usedto discover the rank by monitoring the norms of the remaining unreduced columns andterminating the factorization when the maximum value falls below some relative tolerance.3.5Comparison of MethodsWe have now seen a number of methods for solving least squares problems. The choiceamong them depends on the particular problem being solved and involves trade-offs amongefficiency, accuracy, and robustness.The normal equations method is easy to implement: it simply requires matrix multiplication and Cholesky factorization.

Moreover, reducing the problem to an n × n system isvery attractive when m n. By taking advantage of its symmetry, the formation of thenormal equations matrix AT A requires about n2 m/2 multiplications and a similar numberof additions. Solving the resulting linear system by Cholesky factorization requires aboutn3 /6 multiplications and a similar number of additions. Unfortunately, the normal equations method produces a solution whose relative error is proportional to [cond(A)]2 , and√the required Cholesky factorization can be expected to break down if cond(A) ≈ 1/ machor worse.For solving dense linear least squares problems, the Householder method is generally the most efficient and accurate of the orthogonalization methods.

It requires aboutn2 m − n3 /3 multiplications and a similar number of additions. It can be shown that theHouseholder method produces a solution whose relative error is proportional to cond(A) +krk2 [cond(A)]2 , which is the best that can be expected since this is the inherent sensitivityof the solution to the least squares problem itself. Moreover, the Householder method canbe expected to break down (in the back-substitution phase) only if cond(A) ≈ 1/mach orworse.For nearly square problems, m ≈ n, the normal equations and Householder methodsrequire about the same amount of work.

But for highly overdetermined problems, m n,the Householder method requires about twice as much work as the normal equations method.On the other hand, the Householder method is more accurate and more broadly applicablethan the normal equations method. These advantages may not be worth the additionalcost, however, when the problem is sufficiently well-conditioned that the normal equationsmethod provides adequate accuracy.

For rank-deficient or nearly rank-deficient problems, ofcourse, the Householder method with column pivoting can produce a useful solution whenthe normal equations method would fail outright.3.6Software for Linear Least SquaresTable 3.1 is a list of appropriate routines for solving linear least squares problems, both thosehaving full rank and those that are rank-deficient. Most of the routines listed are based on104CHAPTER 3. LINEAR LEAST SQUARESQR factorization.

Many packages also include software for the singular value decomposition(SVD), which can be used to solve least squares problems, although at greater computationalexpense. The SVD provides a particularly robust method for determining numerical rankand dealing with possible rank deficiency, as we will see in Section 4.5.Table 3.1: Software for linear least squares problemsSourceFactorSolveRank-deficientFMMsvdsvdIMSLlqrrrlqrsllsqrrKMNsqrlssqrlsssvdcLAPACKsgeqrfsormqr/strtrs sgeqpf/stzrqfLawson/Hanson [163] hfths1hftiLINPACKsqrdcsqrslsqrstMATLABqr\svdNAGf01axff04anff04jgfNAPACKqroversing/rsolveaNRqrdcmpqrsolvsvdcmp/svbksbNUMALlsqortdec lsqsolsolovrSLATECsqrdcsqrslllsia/sglss/minfitSOL [279]hredlqrvslvmnlnlsaAs published, qrdcmp and qrsolv handle only square matrices, but they are easilymodified to handle rectangular matrices.Conventional software for solving linear least squares problems Ax ≈ b is sometimesimplemented as a single routine, or it may be split into two routines: one for computing afactorization and another for solving the resulting triangular system.

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

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