Gram-Schmidt process with a non-standard inner product and its application to the approximate inverse preconditioning Miro Rozloˇ zn´ ık joint work with Jiˇ r´ ı Kopal, Alicja Smoktunowicz and Miroslav T˚ uma Institute of Computer Science, Czech Academy of Sciences, Prague, Czech Republic Workshop on Structured Preconditioning and Iterative Methods with Applications, Tsinghua Sanya International Mathematics Forum, Sanya, March 24-28, 2014 Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis TSIMF, Sanya, 2014 1 / 26
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 TSIMF, Sanya, 2014 2 / 26
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 TSIMF, Sanya, 2014 3 / 26
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 O ( 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 TSIMF, Sanya, 2014 4 / 26
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 TSIMF, Sanya, 2014 5 / 26
ON THE WAY FROM THE STANDARD TO THE NONSTANDARD INNER PRODUCT ◮ Axel Ruhe. Numerical aspects of. Gram-Schmidt orthogonalization. of vectors. Lin. Alg. and its Appl.,. 52/53 (591-601), 1983. ◮ T. Ericsson, An analysis of orthogonalization in elliptic norms, to appear. ◮ M. Gulliksson: Backward error analysis for the constrained and weighted linear least squares problem when using the weighted QR factorization. SIAM J. Matrix Anal. Appl. 16(2), 675-687 (1995) M. Gulliksson: On the modified GramSchmidt algorithm for weighted and constrained linear least squares problems. BIT Numer. Math. 35(4), 453-468 (1995) ◮ S.J. Thomas, R.V.M. Zahar: Efficient orthogonalization in the M-norm. Congr. Numer. 80, 23-32 (1991) 36. S.J. Thomas, R.V.M. Zahar, : An analysis of orthogonalization in elliptic norms. Congr. Numer. 86, 193-222 (1992) Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis TSIMF, Sanya, 2014 6 / 26
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 TSIMF, Sanya, 2014 7 / 26
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 ) U − 1 � � ∈ R m , n upper triangular A = I : Z = 0 κ ( U ) = κ ( Z ) = κ 1 / 2 ( A ) Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis TSIMF, Sanya, 2014 8 / 26
Inverse factorization and approximate inverse preconditioning ZZ T = AU − 1 U − T A T = A [ A T BA ] − 1 A T BZZ T : orthogonal projector onto R ( BA ) and orthogonal to R ( A ) ZZ T B : orthogonal projector onto R ( A ) and orthogonal to R ( BA ) A = I square and nonsingular: inverse factorization ZZ T = B − 1 Z T ≈ B − 1 Bx = b , approximate inverse ¯ Z ¯ ↓ Z T B ¯ ¯ Zy = ¯ Z T b , x = ¯ Zy , � ¯ Z T B ¯ Z − I � ≤ ? 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 TSIMF, Sanya, 2014 9 / 26
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 , − u ji z j , j = 1 , . . . , i − 1 i i i z i = z ( i − 1 ) u ii = � z ( i − 1 ) / u ii , � B i i modified Gram-Schmidt ≡ SAINV: u ji = � z ( j − 1 ) , z j � B i classical Gram-Schmidt: u ji = � a i , z j � B AINV algorithm: u ji = � z ( j − 1 ) , a j / u jj � B i Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis TSIMF, Sanya, 2014 10 / 26
CLASSICAL AND MODIFIED GRAM-SCHMIDT ALGORITHMS classical (CGS) modified (MGS) AINV algorithm for i = 1 , . . . , n for i = 1 , . . . , n for i = 1 , . . . , n z ( 0 ) z ( 0 ) z ( 0 ) = a i = a i = a i i i i for j = 1 , . . . , i − 1 for j = 1 , . . . , i − 1 for j = 1 , . . . , i − 1 z ( j ) = z ( j − 1 ) z ( j ) = z ( j − 1 ) − � z ( j − 1 ) z ( j ) = z ( j − 1 ) − � z ( j − 1 ) a j − � a i , z j � B z j , z j � B z j , � B � B i i i i i i i i � z ( j − 1 ) j end end end z i = z ( i − 1 ) / � z ( i − 1 ) z i = z ( i − 1 ) / � z ( i − 1 ) z i = z ( i − 1 ) / � z ( i − 1 ) � B � B � B i i i i i i end end end Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis TSIMF, Sanya, 2014 11 / 26
GRAM-SCHMIDT ALGORITHMS WITH COMPLETE REORTHOGONALIZATION classical (CGS2) modified (MGS2) for i = 1 , . . . , n for i = 1 , . . . , n z ( 0 ) z ( 0 ) = a i = a i i i for k = 1 , 2 for k = 1 , 2 for j = 1 , . . . , i − 1 for j = 1 , . . . , i − 1 z ( j ) = z ( j − 1 ) z ( j ) = z ( j − 1 ) − � z ( j − 1 ) − � a i , z j � B z j , z j � B z j i i i i i end end end end z i = z ( i − 1 ) / � z ( i − 1 ) z i = z ( i − 1 ) / � z ( i − 1 ) � B � B i i i i end end Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis TSIMF, Sanya, 2014 12 / 26
CLASSICAL GRAM-SCHMIDT PROVIDES A CHOLESKY FACTOR Exact arithmetic: a i , a j − � j − 1 � � k = 1 u k , j z k � � u j , i = a i , z j = B u j , j B � a i , a j � B − � j − 1 k = 1 u k , i u k , j = u j , j A T BA = U T U � A − ¯ Z ¯ U � ≤ O ( u ) � ¯ Z �� ¯ U � � ¯ U − ¯ Z T BA � ≤ O ( u ) � A �� B �� ¯ Z � U T ¯ Z T BA ) T ¯ A T BA = ¯ U − (¯ U − ¯ U + A T B ( A − ¯ Z ¯ U ) Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis TSIMF, Sanya, 2014 13 / 26
CLASSICAL GRAM-SCHMIDT PROCESS: THE LOSS OF B -ORTHOGONALITY THE LOSS OF ORTHOGONALITY U T ¯ Z T BA ) T ¯ A T BA = ¯ U − (¯ U − ¯ U + A T B ( A − ¯ Z ¯ U ) U T ¯ U T ¯ A T BA = ¯ Z T B ¯ Z ¯ U + ( A − ¯ Z ¯ U ) T B ¯ Z ¯ U + ¯ Z T B ( A − ¯ Z ¯ U ) + ( A − ¯ Z ¯ U ) T B ( A − ¯ Z ¯ U ) Z T BA ) T ¯ U T ( I − ¯ ¯ Z T B ¯ Z )¯ U = (¯ U − ¯ U − ( A − ¯ Z ¯ U ) T B ¯ Z ¯ U assuming 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 ) Miro Rozloˇ zn´ ık (Prague) Gram-Schmidt orthogonalization: rounding error analysis TSIMF, Sanya, 2014 14 / 26
GRAM-SCHMIDT PROCESS WITH REORTHOGONALIZATION ( B -CGS2) z ( 1 ) = a j − Z j − 1 u ( 1 ) 1 : j − 1 , j = ( I − Z j − 1 Z T j − 1 B ) a j , j z ( 2 ) = z ( 1 ) − Z j − 1 u ( 2 ) j − 1 B ) 2 a j = z ( 1 ) 1 : j − 1 , j = ( I − Z j − 1 Z T j j j � � z ( 2 ) ¯ ¯ u 1 : j − 1 , j � ¯ � � � I − ¯ j − 1 B ¯ Z T Z T Z j − 1 � 2 � j u j , j � j − 1 B ¯ ¯ u j , j 1 / u j , j ≤ σ − 1 min ( U ) = σ − 1 min ( B 1 / 2 A ) , � u 1 : j − 1 , j � / u j , j ≤ κ ( B 1 / 2 A ) , finite precision arithmetic: 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 TSIMF, Sanya, 2014 15 / 26
Recommend
More recommend
