Deturck, Wilf - Lecture Notes on Numerical Analysis (523142), страница 22
Текст из файла (страница 22)
Finally, whenever a pair of different rowsor columns are interchanged we reverse the sign of det. Then det holds the determinantof the input matrix when the forward solution phase has ended.Now we have described the basic modules out of which a general purpose program forlinear equations can be constructed. In the next section we are going to discuss the vexingquestion of roundoff error and how to set the tolerance level below which entries are declaredto be zero. A complete formal algorithm that ties together all of these modules, with controlof rounding error, is given at the end of the next section.Exercises 3.31. Make a test problem for the major program that you’re writing by tracing through asolution the way the computer would:Take one of the systems that appears at the end of section 3.2. Transform it toreduced row echelon form step by step, being sure to carry out a complete search ofthe Southeast rectangle each time, and to interchange rows and columns to bring thelargest element found into the pivot position.
Record the column interchanges in τ ,as described above. Record the status of the matrix C after each major loop so you’llbe able to test your program thoroughly and easily.2. Repeat problem 1 on a system of each major type: inconsistent, unique solution, manysolutions.3. Construct a formal algorithm that will invert a matrix, using no more array spacethan the matrix itself.
The idea is that the input matrix is transformed, a column ata time, into the identity matrix, and the identity matrix is transformed, a column ata time, into the inverse. Why store all of the extra columns of the identity matrix?(Good luck!)3.4 How big is zero?974. Show that a matrix A is of rank one if and only if its entries are of the form Aij = fi gjfor all I and j.−→ := c ∗ −−→ + −−→ applied to a matrix A has the5. Show that the operation −row(i)row(j)row(i)same effect as first applying that same operation to the identity matrix I to get acertain matrix E, and then computing EA.−→6. Show that the operation of scaling −row(i):aik :=aikaiik = 1, .
. . , n(3.3.2)has the same effect as first dividing the ith row of the identity matrix by aii to get acertain matrix E, and then computing EA.7. Suppose we do a complete forward solution without ever searching or interchangingrows or columns. Show that the forward solution amounts to discovering a lowertriangular matrix L and an upper triangular matrix U such that LA = U (think of Las a product of several matrices E such as you found in the preceding two problems).3.4How big is zero?The story of the linear algebra subroutine has just two pieces untold: the first concernshow small we will allow a number to be without calling it zero, and the second concernsthe rearrangement of the output to compensate for interchanges of rows and columns thatare done during the row-echelon reduction.The main reduction loop begins with a search of the rectangle that lies Southeast of thepivot position [i, i], in order to locate the largest element that lives there and to use it forthe next pivot.
If that element is zero, the forward solution halts because the remainingpivot candidates are all zero.But “how zero” do they have to be? Certainly it would be to much to insist, whenworking with sixteen decimal digits, that a number should be exactly equal to zero. Alittle more natural would be to declare that any number that is no larger than the size ofthe accumulated roundoff error in the calculation should be declared to be zero, since ourmicroscope lens would then be too clouded to tell the difference.It is important that we should know how large roundoff error is, or might be. Indeed,if we set too small a threshold, then numbers that “really are” zero will slip through, thecalculation will continue after it should have terminated because of unreliability of thecomputed entries, and so forth.
If the threshold is too large, we will declare numbersto be zero that aren’t, and our numerical solution will terminate too quickly because thecomputed matrix elements will be declared to be unreliable when really they are perfectlyOK.The phenomenon of roundoff error occurs because of the finite size of a computer word.If a word consists of d binary digits, then when two d-digit binary numbers are multiplied98Numerical linear algebratogether, the answer that should be 2d bits long gets rounded off to d bits when it is stored.By doing so we incur a rounding error whose size is at most 1 unit in the (d + 1)st place.Then we proceed to add that answer to other numbers with errors in them, and tomultiply, divide, and so forth, some large number of times.
The accumulation of all of thisrounding error can be quite significant in an extended computation, particularly when agood deal of cancellation occurs from subtraction of nearly equal quantities.The question is to determine the level of rounding error that is present, while thecalculation is proceeding. Then, when we arrive at a stage where the numbers of interestare about the same size as the rounding errors that may be present in them, we had betterhalt the calculation.How can we estimate, during the course of a calculation, the size of the accumulatedroundoff error? There are a number of theoretical a priori estimates for this error, but inany given computation these would tend to be overly conservative, and we would usuallyterminate the calculation too soon, thinking that the errors were worse than they actuallywere.We prefer to let the computer estimate the error for us while it’s doing the calculation.True, it will have to do more work, but we would rather have it work a little harder if theresult will be that we get more reliable answers.Here is a proposal for estimating the accumulated rounding error during the progressof a computation.
This method was suggested by Professors Nijenhuis and Wilf. We carryalong an additional matrix of the same size as the matrix C, the one that has the coefficientsand right-hand sides of the equations that we are solving. In this extra matrix we are goingto keep estimates of the roundoff error in each of the elements of the matrix C.In other words, we are going to keep two whole matrices, one of which will contain thecoefficients and the right-hand sides of the equations, and the other of which will containestimates of the roundoff error that is present in the elements of the first one.At any time during the calculation that we want to know how reliable a certain matrixentry is, we’ll need only to look at the corresponding entry of the error matrix to find out.Let’s call this auxiliary matrix R (as in roundoff).
Initially an element Rij might be aslarge as 2−d |Cij | in magnitude, and of either sign. Therefore, to initialize the R matrix wechoose a number uniformly at random in the interval [−|2−d Cij |, |2−d Cij |], and store it inRij for each i and j.
Hence, to begin with, the matrix R is set to randomly chosen valuesin the range in which the actual roundoff errors lie.Then, as the calculation unfolds, we do arithmetic on the matrix C of two kinds. Weeither scale a row by dividing it through by the pivot element, or we pivot a row againstanother row. In each case let’s look at the effect that the operation has on the correspondingroundoff error estimator in the R matrix.In the first case, consider a scaling operation, in which a certain row is divided by thepivot element.
Specifically, suppose we are dividing row i through by the element Cii , and3.4 How big is zero?99let Rii be the corresponding entry of the error matrix. Then, in view of the fact thatCij + RijCij Rij Rii Cij=+−+ terms involving products of two or more errors (3.4.1)Cii + RiiCii CiiCii2we see that the error entries Rij in the row that is being divided through by Cii should becomputed asRijRii CijRij :=−.(3.4.2)CiiCii2In the second case, suppose we are doing a pivoting operation on the kth row. Then foreach column q we do the operation Ckq := Ckq − t ∗ Ciq , where t = Cki .
Now let’s replaceCkq by Ckq + Rkq , replace Ciq by Ciq + Riq and replace t by t + t0 (where t0 = Rki ). Thensubstitute these expressions into the pivot operation above, and keep terms that are of firstorder in the errors (i.e., that do not involve products of two of the errors).Then Ckq + Rkq is replaced byCkq + Rkq − (t + t0 ) ∗ (Ciq + Riq ) = (Ckq − t ∗ Ciq ) + (Rkq − t ∗ Riq − t0 ∗ Ciq )= (new Ckq ) + (new error Rkq ).(3.4.3)It follows that as a result of the pivoting, the error estimator is updated as follows:Rkq := Rkq − Cki ∗ Riq − Rki ∗ Ciq .(3.4.4)Equations (3.4.2) and (3.4.4) completely describe the evolution of the R matrix.
Itbegins life as random roundoff error; it gets modified along with the matrix elements whoseerrors are being estimated, and in return, we are supplied with good error estimates of eachentry while the calculation proceeds.Before each scaling and pivoting sequence we will need to update the R matrix asdescribed above. Then, when we search the now-famous Southeast rectangle for the newpivot element we accept it if it is larger in absolute value that its corresponding roundoffestimator, and otherwise we declare the rectangle to be identically zero and halt the forwardsolution.The R matrix is also used to check the consistency of the input system.
At the end ofthe forward solution all rows of the coefficient matrix from a certain row onwards are filledwith zeros, in the sense that the entries are below the level of their corresponding roundoffestimator. Then the corresponding right-hand side vector entries should also be zero inthe same sense, else as far as the algorithm can tell, the input system was inconsistent.With typical ambiguity of course, this means either that the input system was “really”inconsistent, or just that rounding errors have built up so severely that we cannot decideon consistency, and continuation of the “solution” would be meaningless.Algorithm matalg(C,r,s,m,n,p,opt,eps). The algorithm operates on the matrix C,which is of dimension r × s, where r = max(m, n) + 1.
It solves p systems of m equationsin n unknowns, unless opt= 1, in which case it will calculate the inverse of the matrix inthe first m = n rows and columns of C.100Numerical linear algebramatalg:=proc(C,r,s,m,n,p,opt,eps)local R,i,j,Det,Done,ii,jj,Z,k,psrank;# if opt = 1 that means inverse is expectedif opt=1 then ident(C,r,s,1,n+1,n,1) fi;# initialize random error matrixR:=matrix(r,s,(i,j)->0.000000000001*(rand()-500000000000)*eps*C[i,j]);# set row permutation to the identityfor j from 1 to n do C[r,j]:=j od;# begin forward solutionDet:=1; Done:=false; i:=0;while ((i<m) and not(Done))do# find largest in SE rectangleZ:=searchmat(C,r,s,i+1,i+1,m,n);ii:=Z[1][1]; jj:=Z[1][2];if abs(Z[2])>abs(R[ii,jj]) theni:=i+1;# switch rowsDet:=Det*switchrow(C,r,s,i,ii,i,s);Z:=switchrow(R,r,s,i,ii,i,s);# switch columnsDet:=Det*switchcol(C,r,s,i,jj,1,r);Z:=switchcol(R,r,s,i,jj,1,r);# divide by pivot elementZ:=scaler(C,R,r,s,i,i);Det:=Det*scale(C,r,s,i,i);for k from i+1 to m do# reduce row k against row iZ:=pivotr(C,R,r,s,i,k,i+1);Z:=pivot(C,r,s,i,k,i+1);od;else Done:=true fi;od;psrank:=i;# end forward solution; begin consistency checkif psrank<m thenDet:=0;for j from 1 to p do# check that right hand sides are 0 for i>psrankZ:=searchmat(C,r,s,psrank+1,n+j,m,n+j);if abs(Z[2])>abs(R[Z[1][1],Z[1][2]]) thenprintf("Right hand side %d is inconsistent",j);return;fi;od;fi;# equations are consistent, do back solution3.4 How big is zero?101for j from psrank to 2 by -1 dofor i from 1 to j-1 doZ:=pivotr(C,R,r,s,j,i,psrank+1);Z:=pivot(C,r,s,j,i,psrank+1);C[i,j]:=0; R[i,j]:=0;od;od;# end back solution, insert minus identity in basisif psrank<n then# fill bottom of basis matrix with -IZ:=ident(C,r,s,psrank+1,psrank+1,n-psrank,-1);# fill under right-hand sides with zeroesfor i from psrank+1 to n do for j from n+1 to s do C[i,j]:=0 od od;# fill under R matrix with zeroesfor i from psrank+1 to n do for j from n-psrank to s do R[i,j]:=0 od od;fi;# permute rows prior to outputZ:=scramb(C,r,s,n);# copy row r of C to row r of Rfor j from 1 to n do R[r,j]:=C[r,j] od;Z:=scramb(R,r,s,n);return(Det,psrank,evalm(R));end;If the procedure terminates successfully, it returns a list containing three items: the firstis the determinant (if there is one), the second is the pseudorank of the coefficient matrix,and the third is the matrix of estimated roundoff errors.
![](https://s.studizba.com/z.php?f=/templates/si/i/noavatar.png&w=100&h=100&t=1"){kind=avatar}
