GRAM-SCHMIDT ORTHOGONALIZATION WITH STANDARD AND NON-STANDARD INNER PRODUCT: ROUNDING ERROR ANALYSIS Miro Rozloˇ zn´ ık joint work with Jiˇ r´ ı Kopal, Ahmed Salam, Alicja Smoktunowicz and Miroslav T˚ uma Institute of Computer Science, Czech Academy of Sciences, Prague, Czech Republic SIAM Conference on Applied Linear Algebra, Valencia, Spain, June 18-22, 2012 Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
STANDARD INNER PRODUCT: GRAM-SCHMIDT PROCESS AS QR ORTHOGONALIZATION A = ( a 1 , . . . , a n ) ∈ R m , n , m ≥ rank ( A ) = n orthogonal basis Q of span ( A ) Q = ( q 1 , . . . , q n ) ∈ R m , n , Q T Q = I A = QR , R upper triangular ( A T A = R T R ) Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
CLASSICAL AND MODIFIED GRAM-SCHMIDT ALGORITHMS finite precision arithmetic: Q T ¯ Q T ¯ ¯ q n ) , ¯ Q � = I n , � I − ¯ Q = (¯ q 1 , . . . , ¯ Q � ≤ ? ◮ classical and modified Gram-Schmidt are mathematically equivalent, but they have ”different” numerical properties ◮ classical Gram-Schmidt can be ”quite unstable”, can ”quickly” lose all semblance of orthogonality ◮ Gram-Schmidt with reorthogonalization : ”two-steps are enough” to preserve the orthogonality to working accuracy Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
GRAM-SCHMIDT PROCESS VERSUS ROUNDING ERRORS ◮ modified Gram-Schmidt: assuming O ( u ) κ ( A ) < 1 Q T ¯ � I − ¯ O ( u ) κ ( A ) Q � ≤ 1 −O ( u ) κ ( A ) Bj¨ orck, 1967, Bj¨ orck, Paige, 1992 ◮ classical Gram-Schmidt: assuming c 1 u κ ( A ) < 1 Q T ¯ O ( u ) κ 2 ( A ) � I − ¯ Q � ≤ 1 −O ( u ) κ ( A ) Giraud, van den Eshof, Langou, R, 2005 Barlow, Smoktunowicz, Langou, 2006 ◮ classical or modified Gram-Schmidt with reorthogonalization : assuming O ( u ) κ ( A ) < 1 Q T ¯ � I − ¯ Q � ≤ O ( u ) Giraud, van den Eshof, Langou, R, 2005 Barlow, Smoktunowicz, 2011 Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
CGS x CHOL QR x MGS x CGS2 0 10 −2 10 −4 10 T Q k || −6 loss of orthgonality || I − Q k 10 −8 10 −10 10 −12 10 −14 10 −16 10 1 2 3 4 5 6 7 8 logarithm of the condition number κ (A k ) Stewart, ”Matrix algorithms” book, p. 284, 1998 Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
ORTHOGONALIZATION WITH A NON-STANDARD INNER PRODUCT B ∈ R m , m symmetric positive definite, inner product �· , ·� B A = [ a 1 , . . . , a n ] ∈ R m , n , m ≥ n = rank ( A ) B -orthogonal basis of the range of A : Z = [ z 1 , . . . , z n ] ∈ R m , n , Z T BZ = I A = ZU , U ∈ R n , n upper triangular A T BA = U T U Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
NON-STANDARD ORTHOGONALIZATION AND STANDARD QR B 1 / 2 A = ( B 1 / 2 Z ) U , Z T BZ = ( B 1 / 2 Z ) T ( B 1 / 2 Z ) = I κ ( Z ) ≪ κ 1 / 2 ( B ) κ ( U ) = κ ( B 1 / 2 A ) ≤ κ 1 / 2 ( B ) κ ( A ) finite precision arithmetic: ¯ z n ) , ¯ Z T B ¯ Z � = I , � I − ¯ Z T B ¯ Z = (¯ z 1 , . . . , ¯ Z � ≤ ? U T ¯ U T ¯ ¯ U ≈ A T BA , � A T BA − ¯ U � ≤ ? Z ¯ ¯ U ≈ A , � A − ¯ Z ¯ U � ≤ ? Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
REFERENCE AND GRAM-SCHMIDT IMPLEMENTATIONS B = V Λ V T , Λ 1 / 2 V T A = QU , Z = V Λ − 1 / 2 Q backward stable eigendecomposition + backward stable QR: � ¯ Z T B ¯ Z − I � ≤ O ( u ) � B �� ¯ Z � 2 z ( 0 ) z ( j ) = z ( j − 1 ) = a i , − α ji z j , j = 1 , . . . , i − 1 i i i z i = z ( i − 1 ) α ii = � z ( i − 1 ) /α ii , � B i i modified Gram-Schmidt ≡ SAINV: α ji = � z ( j − 1 ) , z j � B i classical Gram-Schmidt: α ji = � a i , z j � B AINV algorithm: α ji = � z ( j − 1 ) , a j /α jj � B i Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
LOSS OF B -ORTHOGONALITY IN GRAM-SCHMIDT modified Gram-Schmidt: O ( u ) κ ( B ) κ ( B 1 / 2 A ) < 1 O ( u ) � B �� ¯ Z � 2 κ ( B 1 / 2 A ) � I − ¯ Z T B ¯ Z � ≤ 1 − O ( u ) � B �� ¯ Z � 2 κ ( B 1 / 2 A ) classical Gram-Schmidt and AINV algorithm: O ( u ) κ ( B ) κ ( B 1 / 2 A ) κ ( A ) < 1 O ( u ) � B � 1 / 2 � ¯ Z � κ ( B 1 / 2 A ) κ 1 / 2 ( B ) κ ( A ) � I − ¯ Z T B ¯ Z � ≤ 1 − O ( u ) � B � 1 / 2 � ¯ Z � κ ( B 1 / 2 A ) κ 1 / 2 ( B ) κ ( A ) classical Gram-Schmidt with reorthogonalization: O ( u ) κ 1 / 2 ( B ) κ ( B 1 / 2 A ) < 1 � I − ¯ Z T B ¯ Z � ≤ O ( u ) � B �� ¯ Z �� ¯ Z ( 1 ) � Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
THE LOCAL ERRORS IN A NON-STANDARD INNER PRODUCTS general positive definite B : z ( j − 1 ) z ( j − 1 ) z ( j − 1 ) | fl [ � ¯ , ¯ z j � B ] − � ¯ , ¯ z j � B | ≤ O ( u ) � B �� ¯ �� ¯ z j � i i i z j � 2 z j � 2 | 1 − � ¯ B | ≤ O ( u ) � B �� ¯ diagonal positive ( weight matrix ) B : z ( j − 1 ) z ( j − 1 ) z ( j − 1 ) | fl [ � ¯ , ¯ z j � B ] − � ¯ , ¯ z j � B | ≤ O ( u ) � ¯ � B � ¯ z j � B i i i z j � 2 | 1 − � ¯ B | ≤ O ( u ) Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
DIAGONAL CASE IS SIMILAR TO STANDARD CASE modified Gram-Schmidt: O ( u ) κ ( B 1 / 2 A ) < 1 O ( u ) κ ( B 1 / 2 A ) � I − ¯ Z T B ¯ Z � ≤ 1 − O ( u ) κ ( B 1 / 2 A ) classical Gram-Schmidt and AINV algorithm O ( u ) κ 2 ( B 1 / 2 A ) < 1 O ( u ) κ 2 ( B 1 / 2 A ) � I − ¯ Z T B ¯ Z � ≤ 1 − O ( u ) κ 2 ( B 1 / 2 A ) classical Gram-Schmidt with reorthogonalization: O ( u ) κ ( B 1 / 2 A ) < 1 � I − ¯ Z T B ¯ Z � ≤ O ( u ) Gulliksson, Wedin 1992, Gulliksson 1995 Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
NUMERICAL EXPERIMENTS - EXTREMAL CASES 1. κ 1 / 2 ( B ) ≪ κ ( B 1 / 2 A ) 2. κ ( B 1 / 2 A ) ≤ κ 1 / 2 ( B ) 3. B positive diagonal Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
Problem 9 (Hilbert matrix), κ (B)=1.2e5 0 10 −2 10 −4 10 Loss of orthogonality ||I−Z T BZ|| −6 10 MGS CGS CGS2 −8 10 AINV EIG u κ (B) u κ (B) κ (B 1/2 A) −10 10 u κ (B) κ (B 1/2 A) κ (A) −12 10 −14 10 −16 10 0 2 4 6 8 10 12 14 16 10 10 10 10 10 10 10 10 10 condition number (B 1/2 A) Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
Problem 11 (Hilbert matrix) 0 10 −2 10 −4 10 Loss of orthogonality ||I−Z T BZ|| −6 10 −8 10 −10 10 MGS CGS −12 10 CGS2 AINV EIG −14 u κ (B) 10 u κ (B) κ (B 1/2 A) u κ (B) κ (B 1/2 A) κ (A) −16 10 0 2 4 6 8 10 12 14 16 10 10 10 10 10 10 10 10 10 condition number (B) Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
Problem 3 (diagonal matrix) 0 10 −2 10 −4 10 Loss of orthogonality ||I−Z T BZ|| −6 10 −8 10 MGS CGS −10 10 CGS2 AINV EIG −12 u κ (B 1/2 A) 10 u κ 2 (B 1/2 A) −14 10 −16 10 0 2 4 6 8 10 12 14 16 10 10 10 10 10 10 10 10 10 condition number (B) Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
ELEMENTARY SR FACTORIZATION � � 0 I m ∈ R 2 m , 2 m skew-symmetric and orthogonal J = − I m 0 A = [ a 1 , a 2 ] ∈ R 2 m , 2 , non-isotropic with a T 1 Ja 2 � = 0 � � 0 − 1 V = [ v 1 , v 2 ] ∈ R 2 m , 2 , symplectic V J V = V T JV = I 2 , v T 1 Jv 2 = 1 1 0 � � r 11 r 12 ∈ R 2 , 2 upper triangular with r 11 r 22 = a T A = VR , R = 1 Ja 2 0 r 22 Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
GRAM-SCHMIDT WITH SKEW-SYMMETRIC BILINEAR FORM finite precision arithmetic: V J ¯ ¯ v 2 ) , ¯ v T V = (¯ v 1 , ¯ V � = I , | ¯ 1 J ¯ v 2 − 1 | ≤ ? loss of symplecticity: assuming | a T 1 Ja 2 | > O ( u ) � a 1 �� a 2 � O ( u ) � a 1 �� a 2 � | a T 1 Ja 2 | v T | ¯ 1 J ¯ v 2 − 1 | ≤ 1 − O ( u ) � a 1 �� a 2 � | a T 1 Ja 2 | Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
Thank you for your attention!!! Reference: J. Kopal, R, A. Smoktunowicz, and M. T˚ uma: Rounding error analysis of orthogonalization with a non-standard inner product, to appear in BIT. Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis SIAM Conference on Linear Algebra
Recommend
More recommend
