MATH 612 Computational methods for equation solving and function - PowerPoint PPT Presentation

MATH 612 Computational methods for equation solving and function minimization – Week # 6 F .J.S. Spring 2014 – University of Delaware FJS MATH 612 1 / 58

Plan for this week Discuss any problems you couldn’t solve from previous lectures We will cover Lectures 17, part of 16 and move on to 20 The second coding assignment is due Friday. Radio silence has been decreed for this assignment: you can only discuss it with your coding partner. Remember that... ... I’ll keep on updating, correcting, and modifying the slides until the end of each week. FJS MATH 612 2 / 58

MORE ON BACKWARD STABILITY FJS MATH 612 3 / 58

Review Arithmetic For every real number fl ( x ) = x ( 1 + ǫ ) , | ǫ | ≤ ǫ machine . For every pair of floating point numbers and ∗ ∈ { + , − , × , / } , x ∗ y = ( x ∗ y )( 1 + ǫ ) , | ǫ | ≤ ǫ machine Backward stability An algorithm � f : X → Y to approximate a problem f : X → Y is backward stable, when for all x ∈ X , there exists � x ∈ X such that x ) = � f ( � � � f ( x ) , x − x � ≤ C ǫ machine � x � . FJS MATH 612 4 / 58

The dot product We want to show that the dot product in floating point arithmetic is backward stable. Note that algorithms have to be defined in a fully deterministic way. For instance, given two vectors ( x 1 , . . . , x m ) , ( y 1 , . . . , y m ) , our problem is f : X → R (where X = R n × R n ) given by � f (( x , y )) = f ( x , y ) = x j y j . j The algorithm is, for instance � �� � � f ( x , y ) = . . . ( fl ( x 1 ) × fl ( y 1 )) + ( fl ( x 2 ) × fl ( y 2 )) � � + ( fl ( x 3 ) × fl ( y 3 )) . . . + ( fl ( x m ) × fl ( y m )) FJS MATH 612 5 / 58

Vectors with two components Parturbations arise in floating point representations and arithmetic operations � f ( x , y ) = ( fl ( x 1 ) × fl ( y 1 )) + ( fl ( x 2 ) × fl ( y 2 )) � = x 1 ( 1 + ǫ 1 ) y 1 ( 1 + ǫ 2 ) ( 1 + ǫ 3 ) � �� � � �� � due to × due to fl ( x 1 ) � + x 2 ( 1 + ǫ 4 ) y 2 ( 1 + ǫ 5 )( 1 + ǫ 6 ) ( 1 + ǫ 7 ) � �� � due to + = x 1 ( 1 + ǫ 1 )( 1 + ǫ 7 ) y 1 ( 1 + ǫ 2 )( 1 + ǫ 3 ) � �� � � �� � � � x 1 y 1 + x 2 ( 1 + ǫ 4 )( 1 + ǫ 7 ) y 2 ( 1 + ǫ 5 )( 1 + ǫ 6 ) � �� � � �� � � � x 2 y 2 f ( � x , � = y ) FJS MATH 612 6 / 58

Vectors with two components (cont’d) Using that | ǫ j | ≤ ǫ machine , we cand bound | ( 1 + ǫ i )( 1 + ǫ j ) − 1 | ≤ 2 ǫ machine + O ( ǫ 2 machine ) . so | x i | ( 2 ǫ machine + O ( ǫ 2 | � x i − x i | ≤ machine )) | y i | ( 2 ǫ machine + O ( ǫ 2 | � y i − y i | ≤ machine )) We can wrap up the bounds with the formula y ) − ( x , y ) � ∞ ≤ � ( x , y ) � ∞ ( 2 ǫ machine + O ( ǫ 2 � ( � x , � machine )) FJS MATH 612 7 / 58

Let’s see with three-vectors (11 different ǫ ) � � � f ( x , y ) = ( fl ( x 1 ) × fl ( y 1 )) + ( fl ( x 2 ) × fl ( y 2 )) + ( fl ( x 3 ) × fl ( y 3 )) �� = x 1 ( 1 + ǫ 1 ) y 1 ( 1 + ǫ 2 )( 1 + ǫ 3 ) � + x 2 ( 1 + ǫ 4 ) y 2 ( 1 + ǫ 5 )( 1 + ǫ 6 ) ( 1 + ǫ 7 ) � + x 3 ( 1 + ǫ 8 ) y 3 ( 1 + ǫ 9 )( 1 + ǫ 10 ) ( 1 + ǫ 11 ) = x 1 ( 1 + ǫ 1 )( 1 + ǫ 7 ) y 1 ( 1 + ǫ 2 )( 1 + ǫ 3 )( 1 + ǫ 11 ) � �� � � �� � � � x 1 y 1 + . . . and, collecting carefully, | x i | ( 2 ǫ machine + O ( ǫ 2 | � x i − x i | ≤ machine )) | y i | ( 3 ǫ machine + O ( ǫ 2 | � y i − y i | ≤ machine )) FJS MATH 612 8 / 58

UPPER TRIANGULAR SYSTEMS FJS MATH 612 9 / 58

Back substitution Solve the equations backwards m � r i , j x j = b i , i = m , m − 1 , . . . , 1 j = i by computing � � m � x i = b i − r i , j x j / r i , i , i = m , m − 1 , . . . , 1 . j = i + 1 x=zeros(m,1); for i=m:-1:1 x(i)=(b(i)-r(i,i+1:end)*x(i+1:end))/r(i,i); % row times column end FJS MATH 612 10 / 58

With a picture x b In red what is being used when i = 3. In blue the elements of x that have already been computed, but they are being used in this step. Entries in green are inactive. FJS MATH 612 11 / 58

Recursive back substitution In the previous algorithm, we work equation by equation. Now, we modify the right-hand-side in each step and reduce the size of the system by one. The idea is that once x i is computed, its contribution to all the preceding equations is subtracted. Initial steps: b ( 0 ) = b i i x m = b ( 0 ) m / r m , m b ( 1 ) = b ( 0 ) − r i , m x m i = 1 , . . . , m − 1, i i x m − 1 = b ( 1 ) m / r m − 1 , m − 1 b ( 2 ) = b ( 1 ) − r i , m − 1 x m − 1 i = 1 , . . . , m − 2 , i i x m − 2 = ... In this algorithm we access the upper triangular matrix by columns. Back substitution uses it by rows. FJS MATH 612 12 / 58

Recursive back substitution (2) for i=m:-1:1 x(i)=b(i)/r(i,i); b(1:i-1)=b(1:i-1)-r(1:i-1,i)*x(i); end The flop count for both algorithms is ∼ m 2 . Remark. Forward susbtitution for a lower triangular matrix is back substitution for a system where we count equations and unknowns in the reverse order. Therefore, everything you show for back susbtitution is automatically shown for forward substitution. FJS MATH 612 13 / 58

With a picture x b In red what is being used when i = 3. In blue the elements of x that have already been computed. Entries in green are inactive. FJS MATH 612 14 / 58

BKWD STABILITY OF BACK SUBS FJS MATH 612 15 / 58

In what sense? We have an invertible upper triangular system x = R − 1 b Rx = b , and solve it using back susbtitution and floating point operations. The operator is f ( R ) = x = R − 1 b , that is we consider the matrix to be the input and we fix b . The algorithm is the result of solving with back susbtitution and floating point operations. For simplicity , we will assume that the entries of R and b are already floating point numbers, so we will only care about how the floating point operations in back susbtitution affect the result. FJS MATH 612 16 / 58

In what sense? (2) Backward stability means: given R we can find another upper triangular matrix � R such that R ) = � R − 1 b = f ( � � f ( R ) = � x , where � x is the computed solution, and � � R − R � ≤ C ǫ machine � R � . In other words, we can find � R such that � � � R � x = b = Rx R − R � ≤ C ǫ machine � R � , or also ( R + δ R ) � x = b , � δ R � ≤ C ǫ machine � R � . FJS MATH 612 17 / 58

The 2 × 2 case (Exact) (Computed) x 2 = b 2 / r 22 � x 2 = b 2 / r 22 � x 1 = ( b 1 − ( r 12 × � x 1 = ( b 1 − r 12 x 2 ) / r 11 x 2 )) / r 11 Two tricks to handle ǫ ... if | ǫ | ≤ ǫ machine , then 1 | ǫ ′ | ≤ ǫ machine + O ( ǫ 2 1 + ǫ = machine ) , 1 + ǫ ′ if | ǫ 1 | , | ǫ 2 | ≤ ǫ machine , then | ǫ 3 | ≤ 2 ǫ machine + O ( ǫ 2 ( 1 + ǫ 1 )( 1 + ǫ 2 ) = 1 + ǫ 3 , machine ) . FJS MATH 612 18 / 58

The 2 × 2 case (2) (Exact) (Computed) � x 2 = b 2 / r 22 x 2 = b 2 / r 22 � x 1 = ( b 1 − ( r 12 × � x 1 = ( b 1 − r 12 x 2 ) / r 11 x 2 )) / r 11 x 2 = b 2 b 2 1 ) = b 2 � ( 1 + ǫ 1 ) = r 22 r 22 ( 1 + ǫ ′ � r 22 x 1 = ( b 1 − r 12 � x 2 ( 1 + ǫ 2 ))( 1 + ǫ 3 ) � ( 1 + ǫ 4 ) r 11 = b 1 − r 12 ( 1 + ǫ 2 ) � 4 ) = b 1 − � r 12 � x 2 x 2 r 11 ( 1 + ǫ ′ 3 )( 1 + ǫ ′ � r 11 This is basically it. We have � R � x = b and r ij − r ij | = | r ij | ( c ǫ machine + O ( ǫ 2 | � machine )) , c = 1 or 2 . FJS MATH 612 19 / 58

The general result Given an m × m invertible upper triangular matrix R and the system Rx = b , let � x be the solution computed using back substitution with floating point operations (on a computer satisfying the usual hypotheses). Then there exists an upper triangular matrix � R such that � r ij − r ij | = | r ij | ( m ǫ machine + O ( ǫ 2 R � | � x = b , machine )) . In particular back substitution is backward stable. FJS MATH 612 20 / 58

The general result, made more general What we showed (only for the 2 × 2 case, but extendable for larger cases), does not involve the floating point representation of right-hand-side and matrix. In general, our problem would be f ( R , b ) = R − 1 b . The algorithm computes for decreasing values of i � �� � � � � fl ( r i , i + 1 ) × � fl ( r i , m ) × � x i = fl ( b i ) − x i + 1 − . . . − x m / fl ( r i , i ) Backward stability would mean that we can find a triangular matrix � R and a right-hand-side � b such that � � � � R − R � b − b � � x = � R � b , = O ( ǫ machine ) , = O ( ǫ machine ) , � R � � b � where � x is the computed solution. The proof for the 2 × 2 case is quite simple. You should try and do it in full detail. FJS MATH 612 21 / 58

HOUSEHOLDER REVISITED FJS MATH 612 22 / 58

What Householder triangularization produces Let us start with an invertible matrix A . The Householder method delivers a QR decomposition A = QR , with R given as an upper triangular matrix and Q stored as a lower triangular matrix U , whose columns are the vectors originating the Householder reflectors: Q = ( I − 2 u 1 u ∗ 1 )( I − 2 u 2 u ∗ 2 ) . . . ( I − 2 u m u ∗ m ) ,   0  .  .   . v j ∈ C m − j + 1 . u j =  ,    0 v j FJS MATH 612 23 / 58