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

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

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

Describe how your routine wouldchange if you included partial pivoting. Describe how your routine would change if the81system were positive definite and you computed the Cholesky factorization instead of theLU factorization.2.13 The determinant of a triangular matrixis equal to the product of its diagonal entries.Use this fact to develop a routine for computing the determinant of an arbitrary n × n matrix A by using its LU factorization. You mayuse a library routine for Gaussian eliminationwith partial pivoting to obtain the LU factorization, or you may design your own routine.How can you determine the proper sign for thedeterminant? To avoid risk of overflow or underflow, you may wish to consider computingthe logarithm of the determinant instead of theactual value of the determinant.2.14 Write programs implementing matrixmultiplication C = AB, where A is m × nand B is n × k, in two different ways:(a) Compute the mk inner products of rows ofA with columns of B,(b) Form each column of C as a linear combination of columns of A.In BLAS terminology (see Section 2.7.2), thefirst implementation uses sdot, whereas thesecond uses saxpy.

Compare the performanceof these two implementations on your computer. You may need to try fairly large matrices before the differences in performance become significant. Find out as much as you canabout your computer system (e.g., cache sizeand cache management policy), and use thisinformation to explain the results you observe.2.15 Implement Gaussian elimination usingeach of the six different orderings of the triplenested loop and compare their performance onyour computer.

For purposes of this exercise,you may ignore pivoting for numerical stability, but be sure to use test matrices that do notrequire pivoting. You may need to try a fairlylarge system before the differences in performance become significant. Find out as muchas you can about your computer system (e.g.,cache size and cache management policy), anduse this information to explain the results youobserve.82CHAPTER 2. SYSTEMS OF LINEAR EQUATIONS2.16 Both forward- and back-substitutionfor solving triangular linear systems involvenested loops whose two indices can be takenin either order. Implement both forward- andback-substitution using each of the two index orderings (a total of four algorithms), andcompare their performance for triangular testmatrices of various sizes.

You may need to trya fairly large system before the differences inperformance become significant. Is the bestchoice of index orderings the same for bothalgorithms? Find out as much as you canabout your computer system (e.g., cache sizeand cache management policy), and use thisinformation to explain the results you observe.2.17 Consider a horizontal cantilevered beamthat is clamped at one end but free along theremainder of its length. A discrete model ofthe forces on the beam yields a system of linear equations Ax = b, where the n × n matrixA has the banded form9 −410 ··· ···0..

.. −4.6 −41. .. .. 1 −4.6 −41. , 0 ... ... ... ... ...0  . ... ..1 −46 −41  ... ...1 −45 −2 0 ··· ···01 −21the n-vector b is the known load on the bar(including its own weight), and the n-vector xrepresents the resulting deflection of the barthat is to be determined. We will take the barto be uniformly loaded, with bi = 1 for eachcomponent of the load vector.(a) Letting n = 100, solve this linear system using both a standard library routine fordense linear systems and a library routine designed for band (or more general sparse) systems. How do the two routines compare in thetime required to compute the solution? Howwell do the answers obtained agree with eachother?(b) Verify that the matrix A has the UL factorization A = RRT , where R is an uppertriangular matrix of the form2 −210 ···0..

.. 0.1 −21.  .. . ... .... .....0. .... ..1 −21  ... ...1 −2 0······ ···01Letting n = 1000, solve the linear system usingthis factorization (two triangular solves will berequired). Also solve the system in its originalform using a band solver as in part a. Howwell do the answers obtained agree with eachother? Which approach seems more accurate?What is the condition number of A, and whataccuracy does it suggest that you should expect? Try iterative refinement to see if theaccuracy or residual improves for the less accurate method.Chapter 3Linear Least SquaresWhat meaning should we attribute to a system of linear equations Ax = b if the matrix Ais not square? Since a nonsquare matrix cannot have an inverse, the system of equationsmust have either no solution or a nonunique solution.

Nevertheless, it is often useful todefine a unique vector x that satisfies the linear system in an approximate sense. In thischapter we will see how such problems arise and consider methods for solving them.Let A be an m × n matrix. We will be concerned primarily with the most commonlyoccurring case, m > n, which is called overdetermined because there are more equationsthan unknowns. Such a system usually has no exact solution in the usual sense. Later onwe will also briefly consider the underdetermined case, m < n, with fewer equations thanunknowns.3.1Data FittingPerhaps the most common source of overdetermined linear systems is data fitting, especially when the data have some random error associated with them, as do most empiricallaboratory measurements or other observations of nature.

Given m data points (ti , yi ), wewish to find the n-vector x of parameters that gives the “best fit” to the model functionf (t, x), where f : Rn+1 → R. By best fit we meanminxmX(yi − f (ti , x))2 ,i=1which is called a least squares solution because the sum of squares of differences betweenmodel and data is minimized. Such a problem is usually known as regression analysis instatistics. Note that the quantity being minimized is just the square of the Euclidean 2norm. Other norms, such as the 1-norm or ∞-norm, can be used instead, but they are lessconvenient computationally and give different results with different statistical properties.A least squares problem is linear if the function f is linear in the components of theparameter vector x, which means that f is a linear combinationf (t, x) = x1 φ1 (t) + x2 φ2 (t) + · · · + xn φn (t)8384CHAPTER 3.

LINEAR LEAST SQUARESof functions φj that depend only on t.Example 3.1 Data Fitting. Polynomial fitting, withf (t, x) = x1 + x2 t + x3 t2 + · · · + xn tn−1 ,is a linear least squares problem because a polynomial is linear in its coefficients xj , althoughnonlinear in the independent variable t. An example of a nonlinear least squares data-fittingproblem is a sum of exponentialsf (t, x) = x1 ex2 t + · · · + xn−1 exn t .We will consider nonlinear least squares problems in Section 6.4, but in this chapter wewill confine our attention to linear least squares problems.

We will be concerned only withnumerical algorithms for solving least squares problems. For the many important statisticalconsiderations in formulating least squares problems and in interpreting the results, consultany book on regression analysis or multivariate statistics.3.2Linear Least SquaresA linear least squares data-fitting problem can be written in matrix notation asAx ≈ b,where aij = φj (ti ) and bi = yi . For example, in fitting a quadratic polynomial, which hasthree parameters, to the five data points (t1 , y1 ), . . .

, (t5 , y5 ), the matrix A is 5 × 3, and theproblem has the form y11 t1 t21   y2  1 t2 t22  x1 2  Ax =  1 t3 t3  x2 ≈  y3  = b. y4  1 t4 t2  x342y51 t5 t5A matrix A of this particular form, whose columns (or rows) are successive powers of someindependent variable, is called a Vandermonde matrix .In least squares problems we write Ax ≈ b rather than Ax = b because the “equation”is not usually satisfiable exactly.

The approximate nature of least squares solutions shouldnot disturb us, however, because the goal is to smooth out random errors in the dataand capture the underlying trend. The method of least squares was developed by Gaussfor solving problems in astronomy, particularly determining the orbits of celestial bodiessuch as planets and comets. The elliptical orbit of such a body is determined by fiveparameters, so in principle only five observations of its position should be necessary todetermine the complete orbit. Owing to measurement errors, however, an orbit based ononly five observations would be highly unreliable.

Instead, many more observations are3.3. NORMAL EQUATIONS METHOD85taken and a least squares fit performed in order to smooth out the errors and obtain moreaccurate values for the orbital parameters.As we will see, an m × n linear least squares problem Ax ≈ b has a unique solutionprovided that rank(A) = n (i.e., the columns of A are linearly independent). If rank(A) <n, then A is said to be rank-deficient, and the corresponding linear least squares problemdoes not have a unique solution. We will consider the implications of rank deficiency later,but for now we will assume that A has full rank.Example 3.2 Linear Least Squares Data Fitting. We illustrate linear least squaresby fitting a quadratic polynomial to the following five data points:ty−1.01.0−0.50.50.00.00.50.51.02.0The overdetermined 5 × 3 linear system is therefore1 −1.0 1.0  1.0 1 −0.5 0.25  x1 0.5 Ax =  10.0 0.0  x2 ≈  0.0  = b.10.5 0.25x30.5 11.0 1.02.0The solution to this system, which we will see later how to compute, turns out to bex = [ 0.086 0.40 1.4 ]T , which means that the approximating polynomial isp(t) = 0.086 + 0.4t + 1.4t2 .The resulting curve and the original data points are shown in Fig.

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

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