NUMERICS OF THE GRAM-SCHMIDT ORTHOGONALIZATION PROCESS Miro Rozloˇ zn ´ ık Institute of Computer Science, Czech Academy of Sciences, Prague, Czech Republic email: miro@cs.cas.cz joint results with Luc Giraud, Julien Langou and Jasper van den Eshof but also Alicja Smoktunowicz and Jesse Barlow! Seminarium Pier´ scienie, macierze i algorytmy numeryczne, Politechnika Warszawska, Warszawa, November 13, 2009 1
OUTLINE 1. HISTORICAL REMARKS 2. CLASSICAL AND MODIFIED GRAM-SCHMIDT ORTHOG- ONALIZATION 3. GRAM-SCHMIDT IN FINITE PRECISION ARITHMETIC: LOSS OF ORTHOGONALITY IN THE CLASSICAL GRAM- SCHMIDT PROCESS 4. NUMERICAL NONSINGULARITY AND GRAM-SCHMIDT ALGORITHM WITH REORTHOGONALIZATION 2
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 n A = QR , R upper triangular ( A T A = R T R ) 3
HISTORICAL REMARKS I origination of the QR factorization used for orthogonaliza- tion of functions: J. P. Gram : ¨ Uber die Entwicklung reeller Funktionen in Reihen mittelst der Methode der Kleinsten Quadrate. Journal f. r. a. Math., 94: 41-73, 1883. algorithm of the QR decomposition but still in terms of functions: E. Schmidt : Zur Theorie der linearen und nichlinearen Integralgleichungen. I Teil. Entwicklung willk¨ urlichen Funktionen nach system vorgeschriebener. Mathematische Annalen, 63: 433-476, 1907. name of the QR decomposition in the paper on nonsym- metric eigenvalue problem, rumor: the ”Q” in QR was originally an ”O” standing for orthogonal: J.G.F. Francis : The QR transformation, parts I and II. Computer Journal 4:265-271, 332-345, 1961, 1962. 4
HISTORICAL REMARKS II ”modified” Gram-Schmidt (MGS) interpreted as an elimi- nation method using weighted row sums not as an orthog- onalization technique: P.S. Laplace : Theorie Analytique des Probabilit´ es. Courcier, Paris, third edi- tion, 1820. Reprinted in P.S. Laplace. (Evres Compe´ etes. Gauthier-Vilars, Paris, 1878-1912). ”classical” Gram-Schmidt (CGS) algorithm to solve linear systems of infinitely many solutions: E. Schmidt : ¨ Uber die Aufl¨ osung linearen Gleichungen mit unendlich vielen Unbekanten, Rend. Circ. Mat. Palermo. Ser. 1, 25 (1908), pp. 53-77. first application to finite-dimensional set of vectors: G. Kowalewski : Einfuehrung in die Determinantentheorie. Verlag von Veit & Comp., Leipzig, 1909. 5
CLASSICAL AND MODIFIED GRAM-SCHMIDT ALGORITHMS classical (CGS) modified (MGS) for j = 1 , . . . , n for j = 1 , . . . , n u j = a j u j = a j for k = 1 , . . . , j − 1 for k = 1 , . . . , j − 1 u j = u j − ( a j , q k ) q k u j = u j − ( u j , q k ) q k end end q j = u j / � u j � q j = u j / � u j � end end 6
CLASSICAL AND MODIFIED GRAM-SCHMIDT ALGORITHMS finite precision arithmetic: Q T ¯ Q T ¯ ¯ q n ), ¯ Q � = I n , � I − ¯ Q = (¯ q 1 , . . . , ¯ Q � ≤ ? A � = ¯ Q ¯ R , � A − ¯ Q ¯ R � ≤ ? R ? , cond( ¯ ¯ R ) ≤ ? classical and modified Gram-Schmidt are mathematically equiv- alent, but they have ”different” numerical properties classical Gram-Schmidt can be ”quite unstable”, can ”quickly” lose all semblance of orthogonality 7
ILLUSTRATION, EXAMPLE 1 1 1 σ 0 0 L¨ auchli, 1961, Bj¨ orck, 1967: A = 0 σ 0 0 0 σ κ ( A ) = σ − 1 ( n + σ 2 ) 1 / 2 ≈ σ − 1 √ n , σ ≪ 1 � n + σ 2 σ min ( A ) = σ , � A � = assume first that σ 2 ≤ u , so fl(1 + σ 2 ) = 1 if no other rounding errors are made, the matrices computed in CGS and MGS have the following form: 8
ILLUSTRATION, EXAMPLE 1 0 0 1 0 0 σ − 1 2 − 1 σ − 1 2 − 1 √ √ √ √ 6 2 , 1 − 1 1 0 √ √ 0 √ 0 2 6 2 √ 1 2 0 0 √ 0 0 √ 2 3 √ CGS: (¯ q 3 , ¯ q 1 ) = − σ/ 2, (¯ q 3 , ¯ q 2 ) = 1 / 2, √ MGS: (¯ q 3 , ¯ q 1 ) = − σ/ 6, (¯ q 3 , ¯ q 2 ) = 0 complete loss of orthogonality ( ⇐ ⇒ loss of lin. independence, loss of (numerical) rank ): σ 2 ≤ u (CGS), σ ≤ u (MGS) 9
GRAM-SCHMIDT PROCESS VERSUS ROUNDING ERRORS • modified Gram-Schmidt (MGS): assuming ˆ c 1 uκ ( A ) < 1 Q T ¯ ˆ c 2 uκ ( A ) � I − ¯ Q � ≤ 1 − ˆ c 1 uκ ( A ) Bj¨ orck, 1967 , Bj¨ orck, Paige, 1992 • classical Gram-Schmidt (CGS)? c 2 uκ n − 1 ( A ) Q T ¯ ˜ � I − ¯ c 1 uκ n − 1 ( A ) ? Q � ≤ 1 − ˜ Kielbasinski, Schwettlik, 1994 Polish version of the book, 2nd edition 10
TRIANGULAR FACTOR FROM CLASSICAL GRAM-SCHMIDT VS. CHOLESKY FACTOR OF THE CROSS-PRODUCT MATRIX exact arithmetic: � i − 1 a j , a i − k =1 r k,i q k � � r i,j = a j , q i = r i,i � i − 1 = ( a j , a i ) − k =1 r k,i r k,j r i,i The computation of R in the classical Gram-Schmidt is closely related to the left-looking Cholesky factorization of the cross- product matrix A T A = R T R 11
q i ) + ∆ e (1) ¯ r i,j = fl ( a j , ¯ q i ) = ( a j , ¯ i,j � i − 1 a j , fl ( a i − k =1 ¯ q k ¯ r k,i ) + ∆ e (2) + ∆ e (1) = i i,j ¯ r i,i i − 1 r k,i + ∆ e (3) ¯ r i,i ¯ r i,j = � q k ¯ ¯ a j , a i − i k =1 ( a j , ∆ e (2) ) + ∆ e (1) � � + ¯ r i,i i i,j i − 1 r k,j − ∆ e (1) = ( a i , a j ) − � ¯ r k,i [¯ k,j ] k =1 + ( a j , ∆ e (3) � ( a j , ∆ e (2) ) + ∆ e (1) � ) + ¯ r i,i i i i,j 12
CLASSICAL GRAM-SCHMIDT PROCESS: COMPUTED TRIANGULAR FACTOR � i k =1 ¯ r k,i ¯ r k,j = ( a i , a j ) + ∆ e i,j , i < j R T ¯ [ A T A + ∆ E 1 ] i,j = [ ¯ R ] i,j ! � ∆ E 1 � ≤ c 1 u � A � 2 The CGS process is another way how to compute a backward stable Cholesky factor of the cross-product matrix A T A ! 13
CLASSICAL GRAM-SCHMIDT PROCESS: DIAGONAL ELEMENTS u j = ( I − Q j − 1 Q T j − 1 ) a j � u j � = � a j − Q j − 1 ( Q T j − 1 a j ) � = ( � a j � 2 − � Q T j − 1 a j � 2 ) 1 / 2 u j computing q j = j − 1 a j � ) 1 / 2 , ( � a j �−� Q T j − 1 a j � ) 1 / 2 ( � a j � + � Q T � j � a j � 2 = r 2 k,j + ∆ e j,j , � ∆ e j,j � ≤ c 1 u � a j � 2 k =1 ¯ J. Barlow, A. Smoktunowicz, Langou, 2006 14
CLASSICAL GRAM-SCHMIDT PROCESS: THE LOSS OF ORTHOGONALITY R T ¯ A T A + ∆ E 1 = ¯ R , A + ∆ E 2 = ¯ Q ¯ R Q T ¯ R T ( I − ¯ ¯ Q ) ¯ R = − (∆ E 2 ) T A − A T ∆ E 2 − (∆ E 2 ) T ∆ E 2 + ∆ E 1 assuming c 2 uκ ( A ) < 1 c 3 uκ 2 ( A ) Q T ¯ � I − ¯ Q � ≤ 1 − c 2 uκ ( A ) 15
GRAM-SCHMIDT PROCESS VERSUS ROUNDING ERRORS • modified Gram-Schmidt (MGS): assuming ˆ c 1 uκ ( A ) < 1 Q T ¯ ˆ c 2 uκ ( A ) � I − ¯ Q � ≤ 1 − ˆ c 1 uκ ( A ) Bj¨ orck, 1967, Bj¨ orck, Paige, 1992 • classical Gram-Schmidt (CGS): assuming c 2 uκ ( A ) < 1 c 3 uκ 2 ( A ) Q T ¯ � I − ¯ Q � ≤ 1 − c 2 uκ ( A ) ! Giraud, Van den Eshof, Langou, R, 2006 J. Barlow, A. Smoktunowicz, Langou, 2006 16
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 17
GRAM-SCHMIDT ALGORITHMS WITH COMPLETE REORTHOGONALIZATION classical (CGS2) modified (MGS2) for j = 1 , . . . , n for j = 1 , . . . , n u j = a j u j = a j for i = 1 , 2 for i = 1 , 2 for k = 1 , . . . , j − 1 for k = 1 , . . . , j − 1 u j = u j − ( a j , q k ) q k u j = u j − ( u j , q k ) q k end end end end q j = u j / � u j � q j = u j / � u j � end end 18
GRAM-SCHMIDT PROCESS VERSUS ROUNDING ERRORS • Gram-Schmidt (MGS, CGS): assuming c 1 uκ ( A ) < 1 Q � ≤ c 2 uκ 1 , 2 ( A ) Q T ¯ � I − ¯ 1 − c 1 uκ ( A ) Bj¨ orck, 1967, Bj¨ orck, Paige, 1992 Giraud, van den Eshof, Langou, R, 2005 • Gram-Schmidt with reorthogonalization (CGS2, MGS2): assuming c 3 uκ ( A ) < 1 Q T ¯ � I − ¯ Q � ≤ c 4 u Hoffmann, 1988 Giraud, van den Eshof, Langou, R, 2005 19
ROUNDING ERROR ANALYSIS OF REORTHOGONALIZATION STEP u j = ( I − Q j − 1 Q T j − 1 ) a j , v j = ( I − Q j − 1 Q T j − 1 ) 2 a j � u j � = | r j,j | ≥ σ min ( R j ) = σ min ( A j ) ≥ σ min ( A ) � u j � ≤ κ ( A ), � u j � � a j � � v j � = 1 A + ∆ E 2 = ¯ Q ¯ R , � ∆ E 2 � ≤ c 2 u � A � c 1 uκ ( A ) , � ¯ � a j � u j � κ ( A ) c 2 uκ ( A )] − 1 u j � ≤ v j � ≤ [1 − ˜ � ¯ 1 − ˜ � ¯ 20
FLOPS VERSUS PARALLELISM • classical Gram-Schmidt (CGS): mn 2 saxpys • classical Gram-Schmidt with reorthogonalization (CGS2): 2 mn 2 saxpys • Householder orthogonalization: 2( mn 2 − n 3 / 3) saxpys in parallel environment and using BLAS2, CGS2 may be faster than (plain) MGS! Frank, Vuik, 1999, Lehoucq, Salinger, 2001 21
THANK YOU FOR YOUR ATTENTION! 22
Recommend
More recommend
