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

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

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

To preserve thesolution, however, we will need a new type of transformation to achieve triangular form.903.4.1CHAPTER 3. LINEAR LEAST SQUARESTriangular Least Squares ProblemsAs we did with square linear systems, let us consider a least squares problem having anupper triangular matrix. In the overdetermined case, where m > n, such a problem has theform Rbx≈ 1 ,Ob2where R is an n × n upper triangular matrix and where we have partitioned the right-handside vector b similarly. Then we havekrk22 = kb − Axk22 = kb1 − Rxk22 + kb2 k22 .We have no control over the second term, kb2 k22 , in the foregoing sum, but the first termcan be forced to be zero by choosing x to satisfy the triangular systemRx = b1 ,which can be solved for x by back-substitution.

We have therefore found the least squaressolution x and can also conclude that the minimum sum of squares iskrk22 = kb2 k22 .3.4.2Orthogonal TransformationsReducing a matrix to triangular form via Gaussian elimination is not appropriate for solvingleast squares problems, for such a transformation does not preserve the Euclidean normand hence does not preserve the solution to the problem. We now define a type of lineartransformation that does preserve the Euclidean norm.A matrix Q is said to be orthogonal if its columns are orthonormal, i.e., if QT Q = I,the identity matrix. An orthogonal transformation Q preserves the Euclidean norm of anyvector x, sincekQxk22 = (Qx)T Qx = xT QT Qx = xT x = kxk22 .Orthogonal matrices can transform vectors in various ways, such as rotation or reflection;but they do not change the Euclidean length of a vector.

Hence, they preserve the solutionto a linear least squares problem.Orthogonal matrices are of great importance in many areas of numerical computationbecause their norm-preserving property means that they do not amplify error. Thus, forexample, orthogonal transformations can be used to solve square linear systems withoutthe need for pivoting for numerical stability. Unfortunately, orthogonalization methods aresignificantly more expensive computationally than methods based on Gaussian elimination,so their superior numerical properties come at a price that may or may not be worthwhile,depending on context.3.4.3QR FactorizationGiven an m × n matrix A, with m ≥ n, we seek an m × m orthogonal matrix Q such that R,A=QO3.4. ORTHOGONALIZATION METHODS91where R is n × n and upper triangular.

Such a QR factorization transforms the linear leastsquares problem Ax ≈ b into a triangular least squares problem having the same solutionbecause RRxk2 = kQT b −xk2 .kb − Axk2 = kb − QOOAs with Gaussian elimination, we wish to introduce zeros successively into the matrix A,eventually reaching upper triangular form, but do so using orthogonal transformations. Anumber of methods are possible, including• Householder transformations (elementary reflectors)• Givens transformations (plane rotations)• Gram-Schmidt orthogonalizationWe will focus mainly on the use of Householder transformations, the most popular andgenerally the most effective approach in this context; but we will sketch the other twomethods as well.QR factorization has many other uses besides solving least squares problems.

For example, if we partition Q into Q1 , containing the first n columns, and Q2 , containing theremaining m − n columns, then we have RR= [ Q1 Q2 ]= Q1 R.A=QOOThus, if A has full column rank, so that R is nonsingular, then the columns of Q1 form anorthonormal basis for the range space of A; and the columns of Q2 form an orthonormalbasis for its orthogonal complement, which is the same as the null space of AT (i.e., theset of all vectors x such that AT x = o). Such orthonormal bases are useful in eigenvaluecomputations, optimization, and many other problems, as we will see.3.4.4Householder TransformationsA Householder transformation H is a matrix of the formH =I −2vv T,vT vwhere v is a nonzero vector.

From the definition, we see that H = H T = H −1 , so thatH is both orthogonal and symmetric. Given a vector a, we wish to choose the vector v sothat  α100 Ha =  ...  = α  ...  = αe1 .00Using the formula for H, we haveαe1 = Ha = (I − 2vv T2v T a)a=a−v,vT vvT v92CHAPTER 3. LINEAR LEAST SQUARESand hencevT v.2v T aBut the scalar factor is irrelevant in determining v, since it divides out in the formula forH anyway, so we can takev = a − αe1 .v = (a − αe1 )To preserve the norm, we must have α = ±kak2 , and the sign should be chosen to avoidcancellation.

Another potential numerical difficulty is that the computation of kak2 couldincur unnecessary overflow or underflow if the components of a are very large or very small.Dividing a at the outset by its component of largest magnitude avoids this problem. Again,such a scale factor does not change the resulting transformation H.Example 3.4 Householder Transformation. To illustrate the construction just described, we determine a Householder transformation that annihilates all but the first component of the vector 2a = 1.2Following the foregoing recipe, we choose the vector      212αv = a − αe1 = 1 − α 0 = 1 − 0  ,2020where α = ±kak2 = ±3.

Since a1 is positive, we can avoid cancellation by choosing thenegative sign for α. We therefore have    2−35v= 1 −0 = 1.220To confirm that the Householder transformation performs as expected, we compute   25−3Tv a15   Ha = a − 2 T v =  1  − 21 =0,v v30220which shows that the zero pattern of the result is correct and that the norm is preserved.Note that there is no need to form the matrix H explicitly, as the vector v is all we needto apply H to any vector.Using Householder transformations, we can successively introduce zeros column by column below the diagonal of a matrix A to reduce it to upper triangular form. Each Householder transformation must be applied to the remaining unreduced portion of the matrix,but it will not affect any columns already reduced (and hence the zeros are preserved). Inapplying a Householder transformation H to an arbitrary vector x, we note thatHx = (I − 2vT xvv T)x=x−(2)v,vT vvT v3.4.

ORTHOGONALIZATION METHODS93which is substantially cheaper to compute than a general matrix-vector multiplication andrequires only that we know the vector v.The process just described produces a factorization of the formRHn · · · H1 A =,Owhere R is upper triangular. The product of the successive Householder transformationsHn · · · H1 is itself an orthogonal matrix. Thus, if we takeQT = Hn · · · H1 ,or equivalently, Q = H1T · · · HnT ,thenRA=Q.OHence, we have indeed computed the QR factorization of the matrix A, which we cannow use to solve the linear least squares problem. To preserve the solution, however, wemust also transform the right-hand-side vector b by the same sequence of Householdertransformations. We thus solve the equivalent triangular least squares problemRx ≈ QT b.OFor purposes of solving the linear least squares problem, the product Q of the Householdertransformations need not be explicitly formed.

In most software for this problem, R isstored in the upper triangle of the original array containing A, while the vectors v requiredfor forming the individual Householder transformations are stored in the (now zero) lowertriangular portion of A. (Technically, one additional vector of storage is required, since themain diagonals of both Q and R must be stored.) As we have already seen, Householdertransformations are most easily applied in this form anyway (as opposed to explicit matrixvector multiplication), so the vectors v are all that is needed to solve the original leastsquares problem as well as any subsequent problems having the same matrix but differentright-hand-side vectors.

If Q is needed explicitly for some other reason, however, then it canbe computed by multiplying each Householder transformation in sequence times a matrixthat is initially the identity matrix I, but this computation will require additional storage.Example 3.5 Householder QR Factorization. We illustrate Householder QR factorization by using it to solve 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.250.5 11.0 1.02.0The Householder vector v1 for annihilating the subdiagonal entries of the first column of A94CHAPTER 3. LINEAR LEAST SQUARESis   −2.2363.23611  0   1      v1 = 1 −  0  =  1 .1  0   1 101Applying the resulting Householder transformation H1 yields the transformed matrix andright-hand side−1.789−2.2360−1.118 −0.362  0−0.191 −0.405 0.309 −0.655 H1 A =  0 , H1 b =  −0.862  . 00.809 −0.405−0.362 1.13801.3090.345The Householder vector v2 for annihilating the subdiagonal entries of the second column ofH1 A is  000 −0.191   1.581   −1.772     v2 =  0.309  −  0  =  0.309  . 0.809   0   0.809 1.30901.309Applying the resulting Householder transformation H2 yields−1.789−2.2360−1.118 0.632  01.5810 H2 H1 A =  00−0.725  , H2 H1 b =  −1.035  . −0.816  00−0.589 0.404000.047The Householder vector v3 for annihilating the subdiagonal entries of the third column ofH2 H1 A is  000 0   0   0    −  0.935  =  −1.660  .−0.725v3 =    −0.589   0   −0.589 00.0470.047Applying the resulting Householder transformation H3 yields−2.2360−1.118−1.789 0 0.632 1.5810 H3 H2 H1 A = 00.935  , H3 H2 H1 b =  1.336  . 0 0000.026 0000.337We can now solve the upper triangular system Rx = y, where y consists of the first threecomponents of the transformed right-hand side, by back-substitution to obtain0.086x =  0.400  .1.4293.4.

ORTHOGONALIZATION METHODS3.4.595Givens RotationsHouseholder transformations introduce many zeros in a column at once. Although generallygood for efficiency, this approach can be a bit heavy-handed when greater selectivity isneeded in introducing zeros. For this reason, some algorithms use Givens rotations instead,which introduce zeros one at a time.We seek an orthogonal matrix that annihilates a single given component of a vector.One such orthogonal matrix is a plane rotation, often called a Givens rotation in the contextof QR factorization. Given a 2-vector a = [ a1 a2 ]T , we want to choose scalars c and s,which can be interpreted as the cosine and sine of the angle of rotation, such that c sa1α=,−s ca20pwith c2 + s2 = 1, or, equivalently, α = a21 + a22 .

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

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