Heath - Scientific Computing (523150), страница 32
Текст из файла (страница 32)
Show thatall three of these variations are mathematicallyequivalent (though they may differ markedly infinite-precision arithmetic).3.26 Let v be a nonzero n-vector. The hyperplane normal to v is the (n − 1)-dimensionalsubspace of all vectors y such that v T y = o.A reflector is a linear transformation R suchthat Rx = −x if x is a scalar multiple of v,and Rx = x if v T x = o.
Thus, the hyperplaneacts as a mirror: for any vector, its componentwithin the hyperplane is invariant, whereas itscomponent orthogonal to the hyperplane is reversed.(a) Show that R = 2P − I, where P is the orthogonal projector onto the hyperplane normalto v. Draw a picture to illustrate this resultgeometrically.(b) Show that R is symmetric and orthogonal.(c) Show that the Householder transformationH =I −2vv T,vT vis a reflector.(d ) Show that for any two vectors s and t suchthat s 6= t and ksk2 = ktk2 , there is a reflectorR such that Rs = t.(e) Show that any orthogonal matrix Q is aproduct of reflectors.(f ) Illustrate the previous result by expressingthe plane rotationc s,−s cwhere c2 + s2 = 1, as a product of two reflectors.
For some specific angle of rotation, drawa picture to show the mirrors.Computer Problems3.1 For n = 0, 1, . . . , 5, fit a polynomial ofdegree n by least squares to the following data:ty0.01.01.02.72.05.83.06.64.07.55.09.9Make a plot of the original data points alongwith each resulting polynomial curve (you maymake separate graphs for each curve or a single graph containing all of the curves). Whichpolynomial would you say captures the generaltrend of the data better? Obviously, this isa subjective question, and its answer dependson both the nature of the given data (e.g., theuncertainty of the data values) and the purpose of the fit. Explain your assumptions inanswering.3.2 A common problem in surveying is todetermine the altitudes of a series of pointswith respect to some reference point.
Sincethe measurements are subject to error, moreobservations are taken than are strictly necessary to determine the altitudes, and the resulting overdetermined system is solved in theleast squares sense to smooth out errors. Suppose that there are four points whose altitudesx1 , x2 , x3 , x4 are to be determined. In addition to direct measurements of each xi withrespect to the reference point, measurementsare also taken of each point with respect to allof the others. The resulting set of measurements is as follows:x1 = 2.95,x3 = −1.45,x1 − x2 = 1.23,x1 − x4 = 1.61,x2 − x4 = 0.45,x2 = 1.74,x4 = 1.32,x1 − x3 = 4.45,x2 − x3 = 3.21,x3 − x4 = −2.75.Set up the corresponding least squares systemAx ≈ b and use a library routine, or one ofyour own design, to solve it for the best values112of the altitudes.
How do the computed valuescompare with the direct measurements of thesame quantities?3.3 (a) For a series of matrices A of order n,record the execution times for a library routineto compute the LU factorization of A. Usinga linear least squares routine, or one of yourown design, fit a cubic polynomial to the execution times as a function of n.
To obtainreliable results, use a fairly wide range of values for n, say, in increments of 100 from 100up to several hundred, depending on the speedand available memory of the computer you use.You may obtain more accurate timings by averaging several runs for a given matrix size.The resulting cubic polynomial could be usedto predict the execution time for other valuesof n not tried, such as very large values forn. What is the predicted execution time for amatrix of order 10,000?(b) Try to determine the basic execution rate(in floating-point operations per second, orflops) for your computer by timing a knowncomputation, such as matrix multiplication.You can then use this information to determine the complexity of LU factorization, basedon the polynomial fit to the execution times.After converting to floating-point operations,how does the dominant term compare withthe theoretically expected value of 34 n3 (counting both additions and multiplications)? Tryto explain any discrepancy.
If you use a system that provides operation counts automatically, such as MATLAB or some supercomputers,try this same experiment fitting the operationcounts directly.3.4 (a) Solve the following least squaresproblem using any method you like:0.16 0.10 0.26x 0.17 0.11 1 ≈ 0.28 .x22.02 1.293.31(b) Now solve the same least squares problemagain, but this time use the slightly perturbedright-hand side0.27b = 0.25 .3.33CHAPTER 3. LINEAR LEAST SQUARES(c) Compare your results from parts a and b.Can you explain this difference?3.5 A planet follows an elliptical orbit, whichcan be represented in a Cartesian (x, y) coordinate system by the equationay 2 + bxy + cx + dy + e = x2 .(a) Use a library routine, or one of your owndesign, for linear least squares to determinethe orbital parameters a, b, c, d, e, given thefollowing observations of the planet’s position:xyxy1.020.390.560.150.950.320.440.130.870.270.300.120.770.220.160.130.670.180.010.15In addition to printing the values for the orbital parameters, plot the resulting orbit andthe given data points in the (x, y) plane.(b) This least squares problem is nearly rankdeficient.
To see what effect this has on thesolution, perturb the input data slightly byadding to each coordinate of each data point arandom number uniformly distributed on theinterval [−0.005, 0.005] (see Section 13.5) andsolve the least squares problem with the perturbed data.
Compare the new values for theparameters with those previously computed.What effect does this difference have on theplot of the orbit? Can you explain this behavior?(c) Solve the same least squares problem again,for both the original and the perturbed data,this time using a library routine (or one of yourown design) specifically designed to deal withrank deficiency (by using column pivoting, forexample). Such a routine usually includes asan input parameter a tolerance to be used indetermining the numerical rank of the matrix.Experiment with various values for the tolerance, say, 10−k , k = 1, .
. . , 5. What is the resulting rank of the matrix for each value of thetolerance? Compare the behavior of the twosolutions (for the original and the perturbeddata) with each other as the tolerance and theresulting rank change. How well do the resulting orbits fit the data points as the toleranceand rank vary? Which solution would you regard as better: one that fits the data moreCOMPUTER PROBLEMSclosely, or one that is less sensitive to smallperturbations in the data? Why?3.6 To demonstrate the numerical differencebetween the normal equations method and QRfactorization for linear least squares, we need aproblem that is ill-conditioned and also has asmall residual. We can generate such a problem as follows. We will fit a polynomial ofdegree n − 1,pn−1 (t) = x1 + x2 t + x3 t2 + · · · + xn tn−1 ,to m data points (ti , yi ), m > n.
We chooseti = (i − 1)/(m − 1), i = 1, . . . , m, so that thedata points are equally spaced on the interval [0, 1]. We will generate the correspondingvalues yi by first choosing values for the xj ,say, xj = 1, j = 1, . . . , n, and evaluating theresulting polynomial to obtain yi = pn−1 (ti ),i = 1, . . . , m. We could now see whether wecan recover the xj that we used to generatethe yi , but to make it more interesting, we firstrandomly perturb the yi values to simulate thedata error typical of least squares problems.Specifically, we take yi = yi + (2ui − 1) ∗ ,i = 1, . .
. , m, where each ui is a random number uniformly distributed on the interval [0, 1)(see Section 13.5) and is a small positivenumber that determines the maximum perturbation. If you are using the equivalent ofIEEE double precision, reasonable parametersfor this problem are m = 21, n = 12, and = 10−10 .Having generated the data set (ti , yi ) as justoutlined, we will now compare the two methods for computing the least squares solutionto this polynomial data-fitting problem. First,form the system of normal equations for thisproblem and solve it using a library routine forCholesky factorization. Next, solve the leastsquares system using a library routine for QRfactorization. Compare the two resulting solution vectors x.
For which method is the solution more sensitive to the perturbation we introduced into the data? Which method comescloser to recovering the x that we used to generate the data? Does the difference in solutionsaffect our ability to fit the data points (ti , yi )closely by the polynomial? Why?3.7 Use the augmented system method ofSection 3.3.3 to solve the least squares prob-113lem derived in the previous exercise. The augmented system is symmetric but not positivedefinite, so Cholesky factorization is not applicable, but you can use a symmetric indefiniteor LU factorization.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
