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 = [z1, . . . , zn] - m × n matrix having orthogonal columns with respect to the inner product ·, ·A U upper triangular n × n matrix

Z(0) = ZU, ZTAZ = I (Z(0))TAZ(0) = UTU

Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 2 / 20

Condition number of orthogonal and triangular factor

ZTAZ = (A1/2Z)T(A1/2Z) = I = ⇒ A1/2Z is orthogonal with respect to the standard inner product A1/2Z(0) = (A1/2Z)U is a standard QR factorization

κ(Z) ≪ κ1/2(A) κ(U) = κ(A1/2Z(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 ZTAZ = I

ZZT = Z(0)U−1U−T(Z(0))T = Z(0)[(Z(0))TAZ(0)]−1(Z(0))T AZZT :

• rthogonal projector onto R(AZ(0)) and orthogonal to R(Z(0))

ZZTA :

• rthogonal projector onto R(Z(0)) and orthogonal to R(AZ(0))

important case Z(0) square and nonsingular: inverse factorization

ZZT = 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 gives the approximate inverse ¯ Z¯ ZT ≈ A−1 Ax = b ¯ ZTA¯ Zy = ¯ ZTb , where x = ¯ Zy ¯ U approximates the Cholesky factor of (Z(0))TAZ(0)

the loss of orthogonality ¯ ZTA¯ Z − I the factorization error Z(0) − ¯ Z¯ U Cholesky factorization error (Z(0))TAZ(0) − ¯ UT ¯ U

Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 5 / 20

10 10

2

10

4

10

6

10

8

10

10

10

12

10

14

10

16

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 Loss of orthogonality ||I−ZTAZ|| condition number (A) Problem 1a (inverse Hilbert matrix) = 0 = 1e−14 = 1e−12 = 1e−10 = 1e−8 u * ||A||||Z||2||XTX|| 1e−14 * ||A||||Z||2||XTX|| 1e−12 * ||A||||Z||2||XTX|| 1e−10 * ||A||||Z||2||XTX|| 1e−8 * ||A||||Z||2||XTX|| 2 4 6 8 10 12 14 16 18 20 10

−15

10

−10

10

−5

10 10

5

residual norm ||r|| iteration n A A+1e−14||A||*rand(n,n) A+1e−12||A||*rand(n,n) A+1e−10||A||*rand(n,n) A+1e−8||A||*rand(n,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ΛVT
• 2. QR factorization Λ1/2VTZ(0) = QU
• 3. orthogonal-diagonal-orthogonal matrix multiplication Z = VΛ−1/2Q

backward stable eigendecomposition + backward stable QR:

¯ ZTA¯ Z − I ≤ O(u)A¯ Z2

Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 7 / 20

Gram-Schmidt orthogonalization z(j)

i

= z(j−1)

i

− αjizj zi = z(i−1)

i

/αii, αii = zi−1

i

A

modified Gram-Schmidt (MGS) algorithm ≡ SAINV algorithm: αji = z(j−1)

i

, zjA classical Gram-Schmidt (CGS) algorithm: αji = z(0)

i

, zjA AINV algorithm: oblique projections αji = z(j−1)

i

, z(0)

j

A/αjj

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)

i

= z(j−1)

i

− αji¯ zj αji = z(j−1)

i

,¯ zjA z(j)

i ,¯

zjA = (1 − ¯ zj2

A)z(j−1) i

,¯ zjA z(j)

i

= z(j−1)

i

− ¯ αjizj ¯ αji = fl[z(j−1)

i

, zjA] z(j)

i , zjA =

• fl[z(j−1)

i

, zjA] − z(j−1)

i

, zjA

• zj2

A

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)κ(A1/2Z(0)) < 1 I − ¯ ZTA¯ Z ≤ O(u)A¯ Z2κ(A1/2Z(0)) 1 − O(u)A¯ Z2κ(A1/2Z(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)κ(A1/2Z(0))κ(Z(0)) < 1 I − ¯ ZTA¯ Z ≤ O(u)A1/2¯ Zκ(A1/2Z(0))κ1/2(A)κ(Z(0)) 1 − O(u)A1/2¯ Zκ(A1/2Z(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

z(1)

i

= z(0)

i

−

i−1

P

j=1

α(1)

ji zj,

α(1)

ji

= z(0)

i

, zjA z(2)

i

= z(1)

i

−

i−1

P

j=1

α(2)

ji zj,

α(2)

ji

= z(1)

i

, zjA zi = z(2)

i

/αii, αii = z(2)

i

A

O(u)κ1/2(A)κ(A1/2Z(0)) < 1 I − ¯ ZTA¯ Z ≤ O(u)A¯ Z2

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)

i

, zjA] − z(j−1)

i

, zjA| ≤ O(u)Az(j−1)

i

zj |1 − zj2

A| ≤ O(u)Azj2

diagonal (weight matrix) A : |fl[z(j−1)

i

, zjA] − z(j−1)

i

, zjA| ≤ O(u)z(j−1)

i

AzjA |1 − zj2

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)κ(A1/2Z(0)) < 1 I − ¯ ZTA¯ Z ≤ O(u)κ(A1/2Z(0)) 1 − O(u)κ(A1/2Z(0)) CGS and AINV algorithms: O(u)κ2(A1/2Z(0)) < 1 I − ¯ ZTA¯ Z ≤ O(u)κ2(A1/2Z(0)) 1 − O(u)κ2(A1/2Z(0)) CGS with reorthogonalization: O(u)κ(A1/2Z(0)) < 1 I − ¯ ZTA¯ Z ≤ O(u) [standard inner product A = I; MGS: Bj¨

• rck 1967, Bj¨
• rck 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) ≪ κ(A1/2Z(0))
• 2. κ(A1/2Z(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

10 10

2

10

4

10

6

10

8

10

10

10

12

10

14

10

16

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 Loss of orthogonality ||I−ZTAZ|| condition number (A1/2 Z(0)) Problem 9 (Hilbert matrix), κ(A) = 1.2e5

MGS CGS CGS2 AINV EIG u κ(A) u κ(A) κ(A1/2Z(0)) u κ(A) κ(A1/2Z(0)) κ(Z(0)) Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 16 / 20

10 10

2

10

4

10

6

10

8

10

10

10

12

10

14

10

16

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 Loss of orthogonality ||I−ZTAZ|| condition number (A) Problem 11 (Hilbert matrix) MGS CGS CGS2 AINV EIG u κ(A) u κ(A) κ(A1/2Z(0)) u κ(A) κ(A1/2Z(0)) κ(Z0)

Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 17 / 20

10 10

2

10

4

10

6

10

8

10

10

10

12

10

14

10

16

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 Loss of orthogonality ||I−ZTAZ|| condition number (A) Problem 3 (diagonal matrix) MGS CGS CGS2 AINV EIG u κ(A1/2Z(0)) u κ2(A1/2 Z(0))

Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 18 / 20

10 10

2

10

4

10

6

10

8

10

10

10

12

10

14

10

16

10

−16

10

−14

10

−12

10

−10

10

−8

10

−6

10

−4

10

−2

10 Loss of orthogonality ||I−ZTAZ|| condition number (A) Problem 6 (Hilbert matrix) MGS CGS CGS2 AINV EIG u κ(A) u κ(A) κ(A1/2Z(0)) u κ(A) κ(A1/2Z(0)) κ(Z0)

Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 19 / 20

Thank you for your attention!!!

Reference: J. Kopal, R, A. Smoktunowicz, and M. T˚ uma: Rounding error analysis of

• rthogonalization with a non-standard inner product, submitted 2010.

Miroslav Rozloˇ zn´ ık (ICS CAS) Orthogonalization with a non-standard inner product 10th IMACS 20 / 20