ORTHOGONALIZATION WITH A NON-STANDARD INNER PRODUCT WITH THE - PowerPoint PPT Presentation

ORTHOGONALIZATION WITH A NON-STANDARD INNER PRODUCT WITH THE APPLICATION TO PRECONDITIONING Miroslav 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 10th IMACS International Symposium on Iterative Methods in Scientific Computing, Marrakech, Morocco, May 18-21, 2011 Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 1 / 20

Orthogonalization with a non-standard inner product A symmetric positive definite m × m matrix Z ( 0 ) = [ z ( 0 ) 1 , . . . , z ( 0 ) n ] full column rank m × n matrix A -orthogonal basis of span ( Z ( 0 ) ) : Z = [ z 1 , . . . , z n ] - m × n matrix having orthogonal columns with respect to the inner product �· , ·� A U upper triangular n × n matrix Z ( 0 ) = ZU , Z T AZ = I ( Z ( 0 ) ) T AZ ( 0 ) = U T U Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 2 / 20

Condition number of orthogonal and triangular factor Z T AZ = ( A 1 / 2 Z ) T ( A 1 / 2 Z ) = I = ⇒ A 1 / 2 Z is orthogonal with respect to the standard inner product A 1 / 2 Z ( 0 ) = ( A 1 / 2 Z ) U is a standard QR factorization κ ( Z ) ≪ κ 1 / 2 ( A ) κ ( U ) = κ ( A 1 / 2 Z ( 0 ) ) ≤ κ 1 / 2 ( A ) κ ( Z ( 0 ) ) particular case Z ( 0 ) = I : Z = U − 1 upper triangular m × n matrix κ ( U ) = κ ( Z ) Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 3 / 20

Approximation properties of orthogonal factor Z T AZ = I ZZ T = Z ( 0 ) U − 1 U − T ( Z ( 0 ) ) T = Z ( 0 ) [( Z ( 0 ) ) T AZ ( 0 ) ] − 1 ( Z ( 0 ) ) T AZZ T : orthogonal projector onto R ( AZ ( 0 ) ) and orthogonal to R ( Z ( 0 ) ) ZZ T A : orthogonal projector onto R ( Z ( 0 ) ) and orthogonal to R ( AZ ( 0 ) ) important case Z ( 0 ) square and nonsingular: inverse factorization ZZ T = A − 1 Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 4 / 20

Important application: approximate inverse preconditioning Z T ≈ A − 1 ¯ Z gives the approximate inverse ¯ Z ¯ Ax = b Z T A ¯ ¯ Zy = ¯ Z T b , where x = ¯ Zy ¯ U approximates the Cholesky factor of ( Z ( 0 ) ) T AZ ( 0 ) the loss of orthogonality � ¯ Z T A ¯ Z − I � the factorization error � Z ( 0 ) − ¯ Z ¯ U � Cholesky factorization error � ( Z ( 0 ) ) T AZ ( 0 ) − ¯ U T ¯ U � Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 5 / 20

Problem 1a (inverse Hilbert matrix) 0 10 5 10 A A+1e−14||A||*rand(n,n) −2 A+1e−12||A||*rand(n,n) 10 A+1e−10||A||*rand(n,n) A+1e−8||A||*rand(n,n) −4 10 0 10 Loss of orthogonality ||I−Z T AZ|| −6 10 residual norm ||r|| −8 −5 10 10 � = 0 −10 � = 1e−14 10 � = 1e−12 � = 1e−10 � = 1e−8 −10 10 −12 10 u * ||A||||Z|| 2 ||X T X|| 1e−14 * ||A||||Z|| 2 ||X T X|| 1e−12 * ||A||||Z|| 2 ||X T X|| −14 10 1e−10 * ||A||||Z|| 2 ||X T X|| 1e−8 * ||A||||Z|| 2 ||X T X|| −15 10 −16 10 0 2 4 6 8 10 12 14 16 10 10 10 10 10 10 10 10 10 0 2 4 6 8 10 12 14 16 18 20 condition number (A) iteration n Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 6 / 20

Eigenvalue based implementation - EIG 1. spectral decomposition A = V Λ V T 2. QR factorization Λ 1 / 2 V T Z ( 0 ) = QU 3. orthogonal-diagonal-orthogonal matrix multiplication Z = V Λ − 1 / 2 Q backward stable eigendecomposition + backward stable QR: � ¯ Z T A ¯ Z − I � ≤ O ( u ) � A �� ¯ Z � 2 Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 7 / 20

Gram-Schmidt orthogonalization z ( j ) = z ( j − 1 ) − α ji z j i i z i = z ( i − 1 ) α ii = � z i − 1 /α ii , � A i i modified Gram-Schmidt (MGS) algorithm ≡ SAINV algorithm: α ji = � z ( j − 1 ) , z j � A i classical Gram-Schmidt (CGS) algorithm: α ji = � z ( 0 ) , z j � A i AINV algorithm: oblique projections α ji = � z ( j − 1 ) , z ( 0 ) � A /α jj i j Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 8 / 20

Local errors in the (modified) Gram-Schmidt process z ( j ) = z ( j − 1 ) − α ji ¯ z j i i α ji = � z ( j − 1 ) , ¯ z j � A i � z ( j ) A ) � z ( j − 1 ) z j � 2 i , ¯ z j � A = ( 1 − � ¯ , ¯ z j � A i z ( j ) = z ( j − 1 ) − ¯ α ji z j i i α ji = fl [ � z ( j − 1 ) ¯ , z j � A ] i � z ( j ) � fl [ � z ( j − 1 ) , z j � A ] − � z ( j − 1 ) � � z j � 2 i , z j � A = , z j � A A i i Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 9 / 20

Loss of orthogonality in the MGS algorithm: O ( u ) κ ( A ) κ ( A 1 / 2 Z ( 0 ) ) < 1 O ( u ) � A �� ¯ Z � 2 κ ( A 1 / 2 Z ( 0 ) ) � I − ¯ Z T A ¯ Z � ≤ 1 − O ( u ) � A �� ¯ Z � 2 κ ( A 1 / 2 Z ( 0 ) ) Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 10 / 20

Loss of orthogonality in the CGS and AINV algorithms O ( u ) κ ( A ) κ ( A 1 / 2 Z ( 0 ) ) κ ( Z ( 0 ) ) < 1 O ( u ) � A � 1 / 2 � ¯ Z � κ ( A 1 / 2 Z ( 0 ) ) κ 1 / 2 ( A ) κ ( Z ( 0 ) ) � I − ¯ Z T A ¯ Z � ≤ 1 − O ( u ) � A � 1 / 2 � ¯ Z � κ ( A 1 / 2 Z ( 0 ) ) κ 1 / 2 ( A ) κ ( Z ( 0 ) ) Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 11 / 20

Classical Gram-Schmidt (CGS2) with reorthogonalization i − 1 z ( 1 ) = z ( 0 ) α ( 1 ) α ( 1 ) = � z ( 0 ) − P ji z j , , z j � A i i ji i j = 1 i − 1 z ( 2 ) = z ( 1 ) α ( 2 ) α ( 2 ) = � z ( 1 ) − P ji z j , , z j � A i i ji i j = 1 z i = z ( 2 ) α ii = � z ( 2 ) /α ii , � A i i O ( u ) κ 1 / 2 ( A ) κ ( A 1 / 2 Z ( 0 ) ) < 1 � I − ¯ Z T A ¯ Z � ≤ O ( u ) � A �� ¯ Z � 2 Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 12 / 20

Local errors in the inner product and normalization general positive definite A : | fl [ � z ( j − 1 ) , z j � A ] − � z ( j − 1 ) , z j � A | ≤ O ( u ) � A �� z ( j − 1 ) �� z j � i i i | 1 − � z j � 2 A | ≤ O ( u ) � A �� z j � 2 diagonal ( weight matrix ) A : | fl [ � z ( j − 1 ) , z j � A ] − � z ( j − 1 ) , z j � A | ≤ O ( u ) � z ( j − 1 ) � A � z j � A i i i | 1 − � z j � 2 A | ≤ O ( u ) Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 13 / 20

A diagonal similar to orthogonalization with the standard inner product MGS algorithm: O ( u ) κ ( A 1 / 2 Z ( 0 ) ) < 1 O ( u ) κ ( A 1 / 2 Z ( 0 ) ) � I − ¯ Z T A ¯ Z � ≤ 1 − O ( u ) κ ( A 1 / 2 Z ( 0 ) ) CGS and AINV algorithms: O ( u ) κ 2 ( A 1 / 2 Z ( 0 ) ) < 1 O ( u ) κ 2 ( A 1 / 2 Z ( 0 ) ) � I − ¯ Z T A ¯ Z � ≤ 1 − O ( u ) κ 2 ( A 1 / 2 Z ( 0 ) ) CGS with reorthogonalization: O ( u ) κ ( A 1 / 2 Z ( 0 ) ) < 1 � I − ¯ Z T A ¯ Z � ≤ O ( u ) [standard inner product A = I ; MGS: Bj¨ orck 1967, Bj¨ orck and Paige 1992, CGS: Giraud, Langou, R, van den Eshof 2005, Barlow, Langou, Smoktunowicz 2006, CGS2: Giraud, Langou, R, van den Eshof 2005] [weighted least squares problem; MGS: Gulliksson, Wedin 1992, Gulliksson 1995] Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 14 / 20

Numerical experiments - four extremal cases 1. κ 1 / 2 ( A ) ≪ κ ( A 1 / 2 Z ( 0 ) ) 2. κ ( A 1 / 2 Z ( 0 ) ) ≪ κ 1 / 2 ( A ) 3. A positive diagonal 4. Z ( 0 ) = I Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 15 / 20

Problem 9 (Hilbert matrix), κ (A) = 1.2e5 0 10 −2 10 −4 10 Loss of orthogonality ||I−Z T AZ|| −6 10 MGS CGS CGS2 −8 10 AINV EIG u κ (A) u κ (A) κ (A 1/2 Z (0) ) −10 10 u κ (A) κ (A 1/2 Z (0) ) κ (Z (0) ) −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 (A 1/2 Z (0) ) Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 16 / 20

Problem 11 (Hilbert matrix) 0 10 −2 10 −4 10 Loss of orthogonality ||I−Z T AZ|| −6 10 −8 10 −10 10 MGS CGS −12 10 CGS2 AINV EIG −14 u κ (A) 10 u κ (A) κ (A 1/2 Z (0) ) u κ (A) κ (A 1/2 Z (0) ) κ (Z 0 ) −16 10 0 2 4 6 8 10 12 14 16 10 10 10 10 10 10 10 10 10 condition number (A) Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 17 / 20

Problem 3 (diagonal matrix) 0 10 −2 10 −4 10 Loss of orthogonality ||I−Z T AZ|| −6 10 −8 10 MGS CGS −10 10 CGS2 AINV EIG −12 u κ (A 1/2 Z (0) ) 10 u κ 2 (A 1/2 Z (0) ) −14 10 −16 10 0 2 4 6 8 10 12 14 16 10 10 10 10 10 10 10 10 10 condition number (A) Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 18 / 20